home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 / Ham Radio 2000.iso / ham2000 / satellit / stp / stp.cpp < prev    next >
C/C++ Source or Header  |  1994-12-11  |  145KB  |  5,673 lines
  1. /*  
  2.    Simple Tracking Program (STP)
  3.  
  4.  
  5.    All the Satellite/Sun/Moon logic is based on software 
  6.    developed by James R. Miller (G3RUH):
  7.  
  8.       PLAN13 (Version 1.30 February 24, 1993)
  9.       Sun/Moon (private communication)
  10.       ATTHIS (Version 2.10 June 21, 1994)
  11.  
  12.    Utilizes satellite keplerian data in nasa 2 line format or
  13.    Kepler file format
  14.  
  15.    This code written by C. Harper (KD4QIO) (responsible for errors)
  16.  
  17.    Packet:      kd4qio@k4bft.#hsv.al.usa.na
  18.  
  19.    Internet:    harper@huntsville.sparta.com 
  20.                 kd4qio@amsat.org
  21.  
  22.    Snail Mail:  P.O Box 18786
  23.                 Huntsville Al,35804-8786
  24.  
  25.    Compiliers:  Borland Turbo C++ Version 3.1
  26.                 GNU gcc-2.4.5 
  27.  
  28.    Ideas/Software/Help from James R. Miller (G3RUH)
  29.  
  30.    Suggestions/Debugging by Mark Phillips (N2RPZ)
  31. */
  32.  
  33. #define Version "1.8 beta"
  34. #define MSDOS
  35.  
  36. #define BLANK " "
  37.  
  38. /*  
  39.     BLANKS is number of blank spaces placed prior to a line for
  40.     prediction mode.  Allows three-hole-punch to be used
  41.  
  42.                 1234567890
  43. */
  44. #define BLANKS "    "
  45.  
  46. #define BasicYear 1900.0
  47.  
  48. #include <ctype.h>
  49. #include <math.h>
  50. #include <stdio.h>
  51. #include <stdlib.h>
  52. #include <string.h>
  53. #include <time.h>
  54.  
  55. /*  define certain things based on Unix or MS-DOS platform */
  56.  
  57. #ifdef UNIX
  58.  
  59. #define MaxSat 300
  60.  
  61. #include "stp.unx"
  62.  
  63. #define strcmpi(a,b) strcasecmp(a,b)
  64. #define strncmpi(a,b,n) strncasecmp(a,b,n)
  65.  
  66. #else
  67.  
  68. #define MaxSat 50
  69.  
  70. #include "stp.dos"
  71. #include <bios.h>
  72. #include <conio.h>
  73. #include <dos.h>
  74.  
  75. #endif
  76.  
  77. FILE *KeplerDataFile,*ModeDataFile,*NasaDataFile;
  78. FILE *OutputFile,*RevolutionDataFile,*SelectFile;
  79. FILE *ControlFile;
  80. FILE *PlotFile1,*PlotFile2,*PlotFile3;
  81.  
  82. #ifdef UNIX
  83. void ConstantDoppler(int argv, char *argc[]);
  84. #endif
  85.  
  86. void AddQth(float Latitude, float Longitude,float Altitude, char QthName[]);
  87. void AttHis(int Option,int argv, char *argc[]);
  88. void CalcConstants(void);
  89. void Celestial(int argv, char *argc[]);
  90. void ClearTheScreen(void);
  91. void CodeInfo(void);
  92. void DashLine(FILE *WhereToPrint);
  93. void DataPrint(float Currentday);
  94. void DetermineAosLos(int Relative);
  95. void DoCheckSum(char FileName[],int PrintFlag);
  96. void Doppler(float Frequency);
  97. void EclipseSatellitePredict(int argv,char *argc[]);
  98. void EditKepler(void);
  99. void EditNasa(void);
  100. void FootPrintInit(void);
  101. void FootPrint(void);
  102. void Help(void);
  103. void HourMinuteSecond(float time,FILE *WhereToWrite,int DayOfTheWeek);
  104. void InitConstants(void);
  105. void InitAntenna(void);
  106. void InitQth(void);
  107. void InitSatellite(void);
  108. void InitSun(void);
  109. void MatchQth(char WhatToMatch[]);
  110. void MonthDay(float DN);
  111. void MoonAzimuthElevation(int PrintOut);
  112. void MutualRealTime(int argv,char *argc[]);
  113. void MutualViewing(int argv,char *argc[]);
  114. void ObserverSatelliteGeometry(void);
  115. void PercivedSignalStrength(FILE *WhereToWrite);
  116. void PrintHeader(int OptionFlag,float CurrentDay);
  117. void ReadKepler(void);
  118. void ReadKeps(void);
  119. void ReadLocalQth(void);
  120. void ReadMode(void);
  121. void ReadNasa(void);
  122. void ReadOtherQth(void);
  123. void RealTimeDisplay(int argv,char *argc[],int WhatToShow);
  124. void RunCheckSum(int argv, char *argc[]);
  125. void RevolutionCorrection(void);
  126. void SatelliteDump(void);
  127. void SatellitePosition(void);
  128. void SatellitePredict(int argv,char *argc[],int Option);
  129. void SelectSatellites(void);
  130. void SquintRangeRate(void);
  131. void SunAzimuthElevation(int PrintOut);
  132. void SunData(void);
  133. void UTCLOCTime(void);
  134. void WhatMode(void);
  135.  
  136. char *GetStringFromUser(int *m);
  137.  
  138. int CheckSum(FILE *InFile,int line[]);
  139. int DeltTime(float time);
  140. int GetIntNumberFromUser(int *m);
  141. int GridLocator(char *loc,float *Latitude, float *Longitude);
  142. int InputCheck(int minans,int maxans);
  143. int KeyEvent(void);
  144. int RealEndOfLine(char string[]);
  145.  
  146. float AINT(float xyz);
  147. float ATAN2(float x,float y);
  148. float DayNumber(float Y,float M,float D);
  149. float FindAOS(void);
  150. float FindLOS(void);
  151. float FNatn(float Y,float X);
  152. float GetFloatNumberFromUser(int *m);
  153. float LinInt(int npt,float x[],float y[],float xin);
  154. float LOSAngle(float AzimuthToCheck);
  155. float P0to360(float ans);
  156. float TemperatureOfSky(float Frequency);
  157.  
  158. char *MonthNames[] = {"January","February","March","April","May","June",
  159.    "July","August","September","October","November","December"};
  160.  
  161. char *DayNames[] = {"Sunday","Monday","Tuesday","Wednesday","Thursday",
  162.    "Friday","Saturday"};
  163.  
  164. #define MaxKeyWords 12
  165. char *KeyWords[13] = {"satellite:","catalog number:","epoch time:",
  166.                       "element set:","inclination:","ra of node:",
  167.                       "eccentricity:","arg of perigee:",
  168.                       "mean anomaly:","mean motion:","decay rate:",
  169.                       "epoch rev:","checksum:"};
  170. int LenKeyWords[MaxKeyWords+1];
  171.  
  172. typedef struct {
  173.      int SquintOpt,NumModes,SatNum,EleSet;
  174.      float YE,TE,IN,RA,EC,WP,MA,MM,M2,RV,N0,A0,B0,SI,CI,PC,QD,WD,DC;
  175.      float DE,Beacon,Alat,Alon,ax,ay,az,Period;
  176.      char SatsName[40];
  177.                } SatRec;
  178. SatRec SATS[MaxSat];
  179. int NumberOfSats;
  180.  
  181. int NumberOfSatsSelected,WhichSatelliteSelected[MaxSat],UpDown[MaxSat];
  182.  
  183. typedef struct {
  184.       float RiseDay,RiseHour,RiseMinute,RiseSecond;
  185.       float SetDay,SetHour,SetMinute,SetSecond;
  186.                } RSRec;
  187. RSRec RiseSet[MaxSat];
  188.  
  189. #define MaxModes 16
  190. typedef struct {
  191.     int MinPhase,MaxPhase;
  192.     char ModeStr[20];
  193.                } ModeRec;
  194. ModeRec Modes[MaxSat][MaxModes];
  195.  
  196. #define MaxQth 90
  197. typedef struct {
  198.      float LA,LO,HT,Ux,Uy,Uz,Ex,Ey,Ez,Nx,Ny,Nz,Ox,Oy,Oz;
  199.      float V0x,V0y,V0z;
  200.      char QthLoc[40];
  201.                } QthRec;
  202. QthRec Qth[MaxQth];
  203. int NumberOfQth;
  204.  
  205. #define MaxObsAng 360
  206. float ObsAzimuth[MaxObsAng],ObsAngle[MaxObsAng];
  207. int NumberOfObsAngles;
  208.  
  209. int TRUE;
  210.  
  211. float ONE,ZERO,PI,TWOPI,PIOVR2,D2R,THIRD;
  212. float RE,GM,J2,YM,YT,YG,G0,MAS0,MASD,INS,EQC1,EQC2,EQC3,WW,WE,W0;
  213. float CNS,SNS,FL,RP,ZZ,RE2,AstroUnit;
  214.  
  215. float Year,Month,Day,Hour,Minute,Second;
  216.  
  217. float Sx,Sy,Sz,Vx,Vy,Vz,Ax,Ay,Az;
  218. float Rx,Ry,Rz,RangeToSatellite,U,E,N,Azimuth,Elevation;
  219. float AltitudeOfSatellite,RangeRate,SquintAngle;
  220. float RevolutionNumber,CurrentMAd,CurrentMA;
  221. float AP,RAAN;
  222. float MeanAnomaly;
  223.  
  224. float CXx,CXy,CXz,CYx,CYy,CYz,CZx,CZy,CZz;
  225. float SATx,SATy,SATz,ANTx,ANTy,ANTz,VELx,VELy,VELz;
  226.  
  227. int MonthOld,DayOld;
  228.  
  229. float HourOld,MinuteOld,SecondOld,RevolutionOld;
  230. float TimeSinceEpoch,DayNow,TimeNow;
  231. float LocalDay,LocalYear;
  232.  
  233. float DayStart;
  234.  
  235. float TempDay,TempHour,TempMin,TempSec;
  236.  
  237. float Alat,Alon;
  238.  
  239. float GHAA,GHAE,MRSE,MASE;
  240. float Illumination,SunAzimuth,SunElevation;
  241. float SUNx,SUNy,SUNz,SinOfSunAngle;
  242. float CurrentSunAz,CurrentSunEl,CurrentMoonAz,CurrentMoonEl;
  243. double DistanceToSun,DistanceToMoon;
  244.  
  245. float LO1,LO2,LO3,LO4,LO5,LO6,LO7,LO8,LO9,LO10,LO11,LO12;
  246. float DM0,LA1,LA2,LA3,LA4,LA5,LA6,LA7;
  247. float HP0,HP1,HP2,HP3,HP4,HP5,LM0,LMD,MA0,MAD,F0,FD,D0,DD;
  248. float MoonEpoch;
  249.  
  250. int GMTmon,GMTdaym,GMTdayw,GMTdayy,GMTyear,GMThour,GMTmin,GMTsec;
  251. int LOCmon,LOCdaym,LOCdayw,LOCdayy,LOCyear,LOChour,LOCmin,LOCsec;
  252.  
  253. float JulianDate;
  254.  
  255. int WhichObserver,WhichSatellite;
  256. int MonthPointer,DayPointer;
  257.  
  258. int LengthOfLine;
  259.  
  260. char line1[80],line2[80],line3[80],SatName[40];
  261. char FileName[80],CurrentMode[5],ECL[6];
  262.  
  263. int RunningCounter,RunningCounterOld;
  264. float Step,Minutes,Seconds;
  265.  
  266. float FindStep,FindAngleError,FindStepMinimum;
  267.  
  268. float FrequencyModeB,FrequencyModeS;
  269. float BandWidthModeB,BandWidthModeS;
  270. float ReceiveGainModeB,ReceiveGainModeS;
  271. float ReceiveNoiseModeB,ReceiveNoiseModeS;
  272. float AverageNumberOfUsers;
  273.  
  274. #define MaxFootprint 45
  275. float SinOfA[MaxFootprint],CosOfA[MaxFootprint];
  276. float FootprintLat[MaxFootprint],FootprintLon[MaxFootprint];
  277. int NumberOfFootprint;
  278.  
  279. /*  
  280.    flags....
  281.  
  282.     PlusMinus     0 use 0-360 on azimuth output
  283.                   1 use -180 to 180
  284.  
  285.     PlotFlag      0 no plot files generated
  286.                   1 dump data files time,az,el,lat,lon,alt
  287.  
  288.     GroundTrace   0 plot file contains lat,lon,alt
  289.                   1 plot file contains lat,lon,0
  290. */
  291.  
  292. int PlusMinus,PlotFlag,GroundTrace;
  293.  
  294. char PlotTemp[80];
  295. int NumberOfPlots;
  296.  
  297. #define SGN(x) ((x) < 0 ? (-1.0) : (1.0))
  298. #define ABS(x) ((x) < 0 ? (-(x)) : (x))
  299.  
  300. main(int argv, char *argc[])
  301. {
  302.    int Option;
  303.  
  304. /*  initialize and compute constants */
  305.  
  306.    InitConstants();
  307.    CalcConstants();
  308.  
  309. /*  get local qth info */
  310.  
  311.    ReadLocalQth();
  312.  
  313. /*  if no command line arguments dump out info on code */
  314.  
  315.    if (argv <= 1)CodeInfo();
  316.  
  317. /*  read keplerian elements */
  318.  
  319.    ReadKeps();
  320.  
  321. /*  read mode file */
  322.  
  323.    ReadMode();
  324.  
  325. /*  what does the user want */
  326.  
  327.    switch (*argc[1]){
  328.       case 'A':
  329.          Option = 0;
  330.          DetermineAosLos(Option);
  331.          break;
  332.  
  333.       case 'a':
  334.          Option = 1;
  335.          DetermineAosLos(Option);
  336.          break;
  337.  
  338.       case 'c':
  339.          Celestial(argv,argc);
  340.          break;
  341.  
  342.       case 'C':
  343.          RunCheckSum(argv,argc);
  344.          break;
  345.  
  346.       case 'd':
  347.          SatelliteDump();
  348.          break;
  349.          
  350.       case 'e':
  351.          if (KeplerDataFile == NULL){
  352.             EditNasa();
  353.          }else{
  354.             EditKepler();
  355.          }
  356.          break;
  357.  
  358.       case 'E':
  359.          EclipseSatellitePredict(argv,argc);
  360.          break;
  361.  
  362.       case 'h':
  363.          Help();
  364.          break;
  365.  
  366.       case 'j':
  367.          AttHis(0,argv,argc);
  368.          break;
  369.  
  370.       case 'J':
  371.          AttHis(1,argv,argc);
  372.          break;
  373.  
  374.       case 'm':
  375.          MutualViewing(argv,argc);
  376.          break;
  377.  
  378.       case 'M':
  379.          MutualRealTime(argv,argc);
  380.          break;
  381.  
  382.       case 'p':
  383.          Option=0;
  384.          SatellitePredict(argv,argc,Option);
  385.          break;
  386.       
  387.       case 'P':
  388.          Option=1;
  389.          SatellitePredict(argv,argc,Option);
  390.          break;
  391.       
  392.       case 'r':
  393.          Option = 0;
  394.          RealTimeDisplay(argv,argc,Option);
  395.          break;
  396.  
  397.       case 'R':
  398.          Option = 1;
  399.          RealTimeDisplay(argv,argc,Option);
  400.          break;
  401.  
  402.       case 's':
  403.          Option = 2;
  404.          RealTimeDisplay(argv,argc,Option);
  405.          break;
  406.  
  407.       case 'S':
  408.          SelectSatellites();
  409.          break;
  410.  
  411.       default:
  412.          printf (" Cant Match Option %s\n",argc[1]);
  413.          exit(-1);
  414.    }
  415.  
  416.    return 0;
  417. }
  418. void AddQth(float Latitude, float Longitude, float Altitude, char QthName[])
  419. /* 
  420.   add additional Qth information
  421.  
  422.   Latitude   Degrees
  423.   Longitude  Degrees
  424.   Altitude   Meters
  425. */
  426. {
  427.    float D,SL,CL,Rx,Rz,SO,CO;
  428.  
  429.    NumberOfQth++;
  430.  
  431.    Qth[NumberOfQth].LA = Latitude/D2R;
  432.    Qth[NumberOfQth].LO = Longitude/D2R;
  433.    Qth[NumberOfQth].HT = Altitude/1000.0;
  434.  
  435.    strcpy(Qth[NumberOfQth].QthLoc,QthName);
  436.  
  437.    CL = cos(Qth[NumberOfQth].LA);
  438.    SL = sin(Qth[NumberOfQth].LA);
  439.    CO = cos(Qth[NumberOfQth].LO);
  440.    SO = sin(Qth[NumberOfQth].LO);
  441.  
  442.    D = sqrt(RE2*CL*CL + ZZ*SL*SL);
  443.    Rx = RE2/D + Qth[NumberOfQth].HT;
  444.    Rz = ZZ/D + Qth[NumberOfQth].HT;
  445.  
  446. /* observers unit vectors UP EAST and NORTH in geocentric coords */
  447.  
  448.    Qth[NumberOfQth].Ux = CL*CO;
  449.    Qth[NumberOfQth].Ex = -SO;
  450.    Qth[NumberOfQth].Nx = -SL*CO;
  451.    Qth[NumberOfQth].Uy = CL*SO;
  452.    Qth[NumberOfQth].Ey = CO;
  453.    Qth[NumberOfQth].Ny = -SL*SO;
  454.    Qth[NumberOfQth].Uz = SL;
  455.    Qth[NumberOfQth].Ez = ZERO;
  456.    Qth[NumberOfQth].Nz = CL;
  457.  
  458. /* observers XYZ coords at earths surface */
  459.  
  460.    Qth[NumberOfQth].Ox = Rx*Qth[NumberOfQth].Ux;
  461.    Qth[NumberOfQth].Oy = Rx*Qth[NumberOfQth].Uy;
  462.    Qth[NumberOfQth].Oz = Rz*Qth[NumberOfQth].Uz;
  463.  
  464. /* observers velocity in geocentric coords (vz=0) */
  465.  
  466.    Qth[NumberOfQth].V0x = -Qth[NumberOfQth].Oy*W0;
  467.    Qth[NumberOfQth].V0y =  Qth[NumberOfQth].Ox*W0;
  468.    Qth[NumberOfQth].V0z = ZERO;
  469. }
  470. float AINT(float xyz)
  471. {
  472. /* mimic the aint function */
  473.  
  474.    return xyz-fmod(xyz,ONE);
  475. }
  476. float ATAN2(float x,float y)
  477. /* 
  478.  arc tan function with user controlled option to return +- pi or 0-2pi
  479. */
  480. {
  481.    float temp;
  482.  
  483.    if (PlusMinus){
  484.       temp = atan2(x,y);
  485.    }else{
  486.       temp = atan2(x,y);
  487.       if (temp < ZERO){
  488.          temp = TWOPI + temp;
  489.       }
  490.    }
  491.    return temp;
  492. }
  493. void AttHis(int Option,int argv, char *argc[])
  494. /*  compute current attitude given an inital set */
  495. {
  496.    FILE *AttHisFile;
  497.  
  498.    char dummy[2];
  499.    int i,j,nc;
  500.    float T,W,Q,CW,SW,SI,CI,CQ,SQ,A[4][4],U[4],V[4];
  501.    float SVx,SVy,SVz;
  502.    float StartMonth,StartDay,EndDay,StartYear;
  503.    float TotalDays,C,S,MAS,TAS;
  504.  
  505. /*  
  506.   if user did not include satellite name on command line
  507.   then show what we got and ask...
  508. */
  509.  
  510.    if (argv <= 2){
  511.       for (i=0;i<=NumberOfSats;i++){
  512.          printf (" %s\n",SATS[i].SatsName);
  513.       }
  514.  
  515.       printf (" Which Satellite from the list ?  ");
  516.  
  517.       gets(SatName);
  518.    } else {
  519.       strcpy(SatName,argc[2]);
  520.    }
  521.  
  522.    nc=strlen(SatName);
  523.  
  524. /*  can we match it */
  525.  
  526.    WhichSatellite=-1;
  527.  
  528.    for (i=0;i<=NumberOfSats;i++){
  529.       if (strnicmp(SatName,SATS[i].SatsName,nc) == 0){
  530.          if (WhichSatellite == -1){
  531.             WhichSatellite = i;
  532.          }else{
  533.             printf (" Multiple Names.  Please be more specific \n");
  534.             exit(-1);
  535.          }
  536.       }
  537.    }
  538.  
  539.    if (WhichSatellite < 0){
  540.       printf (" Could Not Match Satellite %s in List\n",SatName);
  541.       exit (-1);
  542.    }else{
  543.       printf (" Matched your input as %s\n",SATS[WhichSatellite].SatsName);
  544.    }
  545.  
  546. /*  ask user for inital values */
  547.    printf (" Initial Attitude <Alon Alat> (AO-10   240 -5): ");
  548.    scanf("%f %f",&Alon,&Alat);
  549.  
  550.    printf (" Date above data valid for (AO-10  1 24 1994): ");
  551.    scanf("%f %f %f",&StartMonth,&StartDay,&StartYear);
  552.  
  553. /*  eat carriage return */
  554.    gets(dummy);
  555.  
  556.    Alon /= D2R;
  557.    Alat /= D2R;
  558.  
  559. /*  initialize sun */
  560.    InitSun();
  561.  
  562.    DayNow = DayNumber(StartYear,StartMonth,StartDay);
  563.    
  564. /* 
  565.   Compute the Plane to Interial Transformation Matrix A[][] 
  566.  
  567.   use data based on current satellite epoch
  568. */
  569.  
  570.    T = (DayNow-SATS[WhichSatellite].DE)-SATS[WhichSatellite].TE;
  571.    W = SATS[WhichSatellite].WP + SATS[WhichSatellite].WD*T;
  572.    Q = SATS[WhichSatellite].RA + SATS[WhichSatellite].QD*T;
  573.  
  574.    CW = cos(W);
  575.    SW = sin(W);
  576.    SI = sin(SATS[WhichSatellite].IN);
  577.    CI = cos(SATS[WhichSatellite].IN);
  578.    CQ = cos(Q);
  579.    SQ = sin(Q);
  580.  
  581.    W  = FNatn(SW,CW);
  582.    Q  = FNatn(SQ,CQ);
  583.  
  584.    A[1][1] = CW*CQ - SW*CI*SQ;
  585.    A[1][2] =-SW*CQ - CW*CI*SQ;
  586.    A[1][3] = SI*SQ;
  587.    A[2][1] = CW*SQ + SW*CI*CQ;
  588.    A[2][2] =-SW*SQ + CW*CI*CQ;
  589.    A[2][3] =-SI*CQ;
  590.    A[3][1] = SW*SI;
  591.    A[3][2] = CW*SI;
  592.    A[3][3] = CI;
  593.  
  594. /* Determine attitude in Celestial coordinates AttPTOC()*/
  595.  
  596.    U[1] = cos(Alat)*cos(Alon);
  597.    U[2] = cos(Alat)*sin(Alon);
  598.    U[3] = sin(Alat);
  599.  
  600. /* Plane to Celestial Transform PTOC()*/
  601.  
  602.    for (i=1;i<=3;i++){
  603.       V[i] = ZERO;
  604.       for (j=1;j<=3;j++){
  605.          V[i] = V[i] + A[i][j]*U[j];
  606.       }
  607.    } 
  608. /* end of PTOC() */
  609.  
  610.    SVx = V[1];
  611.    SVy = V[2];
  612.    SVz = V[3];
  613.  
  614. /* finished with AttPTOC() */
  615.  
  616. /*  
  617.   case 0 is compute for today
  618.   case 1 is compute from today forward
  619. */
  620.  
  621. /*  get current time */
  622.    UTCLOCTime();
  623.  
  624.    Year  = (float)GMTyear + BasicYear;
  625.    Month = (float)GMTmon  + ONE;
  626.    Day   = (float)GMTdaym;
  627.  
  628.    DayNow = DayNumber(Year,Month,Day);
  629.  
  630. /* 
  631.   Compute the Plane to Interial Transformation Matrix A[][] 
  632.  
  633.   use data based on current satellite epoch
  634.  
  635.   compute for today 
  636. */
  637.  
  638.    T = (DayNow-SATS[WhichSatellite].DE)-SATS[WhichSatellite].TE;
  639.    W = SATS[WhichSatellite].WP + SATS[WhichSatellite].WD*T;
  640.    Q = SATS[WhichSatellite].RA + SATS[WhichSatellite].QD*T;
  641.  
  642.    CW = cos(W);
  643.    SW = sin(W);
  644.    SI = sin(SATS[WhichSatellite].IN);
  645.    CI = cos(SATS[WhichSatellite].IN);
  646.    CQ = cos(Q);
  647.    SQ = sin(Q);
  648.  
  649.    W  = FNatn(SW,CW);
  650.    Q  = FNatn(SQ,CQ);
  651.  
  652.    A[1][1] = CW*CQ - SW*CI*SQ;
  653.    A[1][2] =-SW*CQ - CW*CI*SQ;
  654.    A[1][3] = SI*SQ;
  655.    A[2][1] = CW*SQ + SW*CI*CQ;
  656.    A[2][2] =-SW*SQ + CW*CI*CQ;
  657.    A[2][3] =-SI*CQ;
  658.    A[3][1] = SW*SI;
  659.    A[3][2] = CW*SI;
  660.    A[3][3] = CI;
  661.  
  662. /*  Attitude in Celestial to Plane AttCTOP() */
  663.  
  664.    U[1] = SVx;
  665.    U[2] = SVy;
  666.    U[3] = SVz;
  667.  
  668. /*  Celestial to Plane CTOP() */
  669.  
  670.    for (i=1;i<=3;i++){
  671.       V[i] = ZERO;
  672.       for (j=1;j<=3;j++){
  673.          V[i] = V[i] + A[j][i]*U[j];
  674.       }
  675.    } 
  676. /*  end of CTOP() */
  677.  
  678.    Alat = asin(V[3]);
  679.    Alon = FNatn(V[2],V[1]);
  680.  
  681.    switch (Option){
  682.       case 0:
  683.          TimeSinceEpoch = (DayNow  - SATS[WhichSatellite].DE) + 
  684.                           (ZERO    - SATS[WhichSatellite].TE);       
  685.  
  686. /*  ma of sun round its orbit */
  687.          MAS = MASE + (MASD*TimeSinceEpoch)/D2R; 
  688.  
  689. /*  suns true anomaly */
  690.          TAS = MRSE + WW*TimeSinceEpoch
  691.               +EQC1*sin(MAS)
  692.               +EQC2*sin(2.0*MAS)
  693.               +EQC3*sin(3.0*MAS);
  694.  
  695.          C = cos(TAS);
  696.          S = sin(TAS);
  697.  
  698. /*  sun unit vector celestial coords....assuming circular orbit */
  699.          SUNx = C;
  700.          SUNy = S*CNS;
  701.          SUNz = S*SNS;
  702.  
  703. /*  sin of sun angle */
  704.          SinOfSunAngle = -(SVx*SUNx + SVy*SUNy + SVz*SUNz);
  705.  
  706. /*  illumination */
  707.  
  708.          Illumination = 100.0*sqrt(ONE-SinOfSunAngle*SinOfSunAngle);
  709.  
  710.          printf (" Current Attitude: Alat/Alon %6.1f/%6.1f",
  711.             Alat*D2R,
  712.             Alon*D2R);
  713.  
  714.          printf (" Sun Angle = %6.1f Illum = %6.1f\n",
  715.             -asin(SinOfSunAngle)*D2R,
  716.             Illumination);
  717.          break;
  718.  
  719.       case 1:
  720.  
  721.          printf (" How Many Days Ahead to Predict: ");
  722.          TotalDays=GetFloatNumberFromUser(&nc);
  723.  
  724.          printf (" Step Size (Days): ");
  725.          Step=GetFloatNumberFromUser(&nc);
  726.  
  727.          StartDay = DayNow;
  728.          EndDay   = StartDay + TotalDays;
  729.  
  730.          printf (" Filename to write output to <cr> for screen: ");
  731.  
  732.          strcpy(FileName,GetStringFromUser(&nc));
  733.  
  734.          if (nc > 0){
  735.             if ((OutputFile = fopen(FileName,"w")) == 0){
  736.                printf ("Can't open file %s\n",FileName);
  737.                exit(-1);
  738.             }
  739.             fprintf (OutputFile,"        Date         Amsat");
  740.             fprintf (OutputFile,"   Alat   Alon   SunAng  Illu\n");
  741.          }else{
  742.             OutputFile = stdout;
  743.          }
  744.  
  745.          AttHisFile = fopen ("atthis.dat","w");
  746.  
  747.          fprintf (AttHisFile,"5\nAmsat Day\nAlat\nAlon\n");
  748.          fprintf (AttHisFile,"Sun Angle\nIllumination\n");
  749.  
  750.          for (DayNow=StartDay;DayNow<=EndDay;DayNow=DayNow+Step){
  751.  
  752.             T = (DayNow-SATS[WhichSatellite].DE)-SATS[WhichSatellite].TE;
  753.             W = SATS[WhichSatellite].WP + SATS[WhichSatellite].WD*T;
  754.             Q = SATS[WhichSatellite].RA + SATS[WhichSatellite].QD*T;
  755.  
  756.             CW = cos(W);
  757.             SW = sin(W);
  758.             SI = sin(SATS[WhichSatellite].IN);
  759.             CI = cos(SATS[WhichSatellite].IN);
  760.             CQ = cos(Q);
  761.             SQ = sin(Q);
  762.  
  763.             W  = FNatn(SW,CW);
  764.             Q  = FNatn(SQ,CQ);
  765.          
  766.             A[1][1] = CW*CQ - SW*CI*SQ;
  767.             A[1][2] =-SW*CQ - CW*CI*SQ;
  768.             A[1][3] = SI*SQ;
  769.             A[2][1] = CW*SQ + SW*CI*CQ;
  770.             A[2][2] =-SW*SQ + CW*CI*CQ;
  771.             A[2][3] =-SI*CQ;
  772.             A[3][1] = SW*SI;
  773.             A[3][2] = CW*SI;
  774.             A[3][3] = CI;
  775.  
  776. /*  Attitude in Celestial to Plane AttCTOP() */
  777.  
  778.             U[1] = SVx;
  779.             U[2] = SVy;
  780.             U[3] = SVz;
  781.          
  782.          /*  Celestial to Plane CTOP() */
  783.  
  784.             for (i=1;i<=3;i++){
  785.                V[i] = ZERO;
  786.                for (j=1;j<=3;j++){
  787.                   V[i] = V[i] + A[j][i]*U[j];
  788.                }
  789.             } 
  790. /*  end of CTOP() */
  791.  
  792.             Alat = asin(V[3]);
  793.             Alon = FNatn(V[2],V[1]);
  794.  
  795.             MonthDay(DayNow);
  796.  
  797.             TimeSinceEpoch = (DayNow  - SATS[WhichSatellite].DE) + 
  798.                              (ZERO    - SATS[WhichSatellite].TE);       
  799.  
  800. /*  ma of sun round its orbit */
  801.             MAS = MASE + (MASD*TimeSinceEpoch)/D2R; 
  802.  
  803. /*  suns true anomaly */
  804.             TAS = MRSE + WW*TimeSinceEpoch
  805.                  +EQC1*sin(MAS)
  806.                  +EQC2*sin(2.0*MAS)
  807.                  +EQC3*sin(3.0*MAS);
  808.  
  809.             C = cos(TAS);
  810.             S = sin(TAS);
  811.  
  812. /*  sun unit vector celestial coords....assuming circular orbit */
  813.             SUNx = C;
  814.             SUNy = S*CNS;
  815.             SUNz = S*SNS;
  816.  
  817. /*  sin of sun angle */
  818.             SinOfSunAngle = -(SVx*SUNx + SVy*SUNy + SVz*SUNz);
  819.  
  820. /*  illumination */
  821.  
  822.             Illumination = 100.0*sqrt(ONE-SinOfSunAngle*SinOfSunAngle);
  823.  
  824.             fprintf (OutputFile," %9s %2.0f %4.0f",
  825.                MonthNames[MonthPointer],
  826.                LocalDay,
  827.                LocalYear);
  828.  
  829.             fprintf (OutputFile," %8.1f %6.1f %6.1f %6.1f %6.1f\n",
  830.                DayNow - 722100.0,
  831.                Alon*D2R,
  832.                Alat*D2R,
  833.                -asin(SinOfSunAngle)*D2R,
  834.                Illumination);
  835.  
  836.             fprintf (AttHisFile," %8.1f %6.1f %6.1f %6.1f %6.1f\n",
  837.                DayNow - 722100.0,
  838.                Alon*D2R,
  839.                Alat*D2R,
  840.                -asin(SinOfSunAngle)*D2R,
  841.                Illumination);
  842.          }
  843.  
  844.          fclose (OutputFile);
  845.          fclose (AttHisFile);
  846.  
  847.          break;
  848.    }
  849. }
  850. void CalcConstants(void)
  851. {
  852.    WW = TWOPI/YT;       /* earths rotation rate, radians/whole day */
  853.    WE = TWOPI + WW;     /* earths rotation rate, radians/day */
  854.    W0 = WE/86400.0;     /* earths rotation rate, radians/sec */
  855.  
  856.    G0 /=D2R;
  857.  
  858.    INS /=D2R;
  859.    CNS = cos(INS);
  860.    SNS = sin(INS);
  861.  
  862.    FL = ONE/298.257;    /* IAU-76 earths ellipsoid */
  863.    RP = RE*(ONE-FL);
  864.    ZZ = RP*RP;
  865.    RE2= RE*RE;
  866. }
  867. void Celestial(int argv, char *argc[])
  868. /*  return the sun/moon elevation and azimuth */
  869. {
  870.    float TEG,DE,TE;
  871.  
  872.    WhichObserver = 0;
  873.  
  874. /*  initailize Qth parameters */
  875.  
  876.    InitQth();
  877.  
  878.    TEG = -1.0;
  879.  
  880. /*  get update rate */
  881.  
  882.    if (argv == 3){
  883.       Step = atof(argc[2]);
  884.    }else{
  885.       Step = 1.0;
  886.    }
  887.  
  888.    while (bioskey(1) == 0){
  889.       UTCLOCTime();
  890.       
  891.       Year  = (float)GMTyear + BasicYear;
  892.       Month = (float)GMTmon  + ONE;
  893.       Day   = (float)GMTdaym;
  894.  
  895.       DayNow = DayNumber(Year,Month,Day);
  896.  
  897.       Hour   = (float)GMThour;
  898.       Minute = (float)GMTmin;
  899.       Second = (float)GMTsec;
  900.  
  901.       TimeNow = (Hour + Minute/60.0 + Second/3600.0)/24.0;
  902.  
  903. /*  if TEG is undefined then bring Sun and Moon data to current time */
  904.  
  905.       if (TEG < ZERO){
  906.          DE  = DayNow;
  907.          TE  = TimeNow;
  908.          TEG = DE-DayNumber(YG,ONE,ZERO)+TE;
  909.  
  910.          GHAE = G0 + TEG*WE ;
  911.          MRSE = G0 + TEG*WW + PI;
  912.          MASE = (MAS0 + MASD*TEG)/D2R;
  913.       }
  914.  
  915.       TimeSinceEpoch = (DayNow  - DE) + (TimeNow - TE);
  916.  
  917. /*  rotate everything */
  918.  
  919.       GHAA = GHAE + WE*TimeSinceEpoch;
  920.  
  921.       if (DeltTime(Step)){
  922.          ClearTheScreen();
  923.  
  924.          printf (" %02.0f:%02.0f:%02.0f",
  925.            Hour,
  926.            Minute,
  927.            Second);
  928.  
  929. /*  get the sun */
  930.  
  931.          SunAzimuthElevation(1);
  932.          printf ("\n");
  933. /*  
  934.   print blanks equivalent to %02.0f:%02.0f:%02.0f
  935.                   bxx:xx:xx
  936. */
  937.          printf ("         ");
  938.  
  939. /* get the moon */
  940.  
  941.          MoonAzimuthElevation(1);
  942.          printf ("\n");
  943.  
  944. /*  save current time for output determination */
  945.  
  946.          HourOld   = Hour;
  947.          MinuteOld = Minute;
  948.          SecondOld = Second;
  949.       }
  950.    }
  951. }
  952. int CheckSum(FILE *InFile,int line[])
  953. /*  determine checksum */
  954. {
  955.    int nlen,idum,ic,i;
  956.  
  957.    idum = 0;
  958.    nlen = 0;
  959.  
  960. /*  
  961.   fgetc returns the ascii code for each character it reads in
  962.  
  963.   68 characters in a line (reads up to checksum value)
  964. */
  965.  
  966.    for (i=1;i<=68;i++){
  967.       ic = fgetc(InFile);
  968.  
  969.       if (ic == -1) {
  970.          return 0;
  971.       }
  972.  
  973.       line[nlen]=ic;
  974.       nlen++;
  975.  
  976. /*  
  977.   if a number....48 == 0   57 == 9 then convert to the numeric value
  978. */
  979.       if (ic >= 48 && ic <= 57){
  980.          idum=idum+ic-48;
  981.       }
  982.  
  983. /*  
  984.    if a + then add 2
  985.    if a - then add 1
  986. */
  987.       if (ic == 43)idum += 2;
  988.       if (ic == 45)idum += 1;
  989.    }
  990.  
  991. /*  mod 10 */
  992.    idum=idum-10*(idum/10);
  993.  
  994. /*  read checksum value and convert from ascii to number */
  995.  
  996.    ic = fgetc(InFile);
  997.  
  998.    line[nlen]=ic;
  999.    nlen++;
  1000.  
  1001. /*  read last character on line (carriage return) */
  1002.    ic=fgetc(InFile);
  1003.  
  1004.    line[nlen]=ic;
  1005.  
  1006.    return idum;
  1007. }
  1008. void ClearTheScreen(void)
  1009. /*  clear screen utility */
  1010. {
  1011. #ifdef UNIX
  1012.    system ("clear");
  1013. #else
  1014.    clrscr();
  1015. #endif
  1016. }
  1017. void CodeInfo(void)
  1018. /*  program information */
  1019. {
  1020. /*  send bells */
  1021.    printf ("\a\a\a\a\n");
  1022.    ClearTheScreen();
  1023.  
  1024.    printf (" Simple Tracking Program (STP) Version: %s\n",Version);
  1025.    printf ("\n");
  1026.    printf (" Propagation model, squint calculation, Sun");
  1027.    printf (" illumination logic\n");
  1028.    printf (" Sun/Moon position model based on software developed");
  1029.    printf (" by G3RUH.\n");
  1030.    printf ("\n");
  1031.    printf (" This implimentation by KD4QIO.\n");
  1032.    printf ("\n");
  1033.    printf (" STP is run using command line options:\n");
  1034.    printf ("\n");
  1035.    printf ("     STP <option> <parameter>\n");
  1036.    printf ("     Options are CASE sensitive... r and R are different.\n");
  1037.    printf ("\n");
  1038.    printf (" A list of options available are shown next.\n");
  1039.    printf ("\n");
  1040.    
  1041.    printf (" Press ENTER to continue\n");
  1042.    gets (line1);
  1043.    ClearTheScreen();
  1044.  
  1045.    printf ("  NONE Gets this list\n");
  1046.    printf ("\n");
  1047.    printf ("  a....AOS/LOS times (relative time) \n");
  1048.    printf ("  A....AOS/LOS times (absolute time) \n");
  1049.    printf ("  c....Sun/Moon real time display\n");
  1050.    printf ("  C....Perform Check Sum calculations on 2-line nasa file\n");
  1051.    printf ("           (can put filename on command line... STP C nasa.dat)\n");
  1052.    printf ("  d....Dump information about satellites (beacon freq,");
  1053.    printf (" alat,alon, # of modes)\n");
  1054.    printf ("  e....View/Edit orbital parameters\n");
  1055.    printf ("  h....Simple suggestions\n");
  1056.    printf ("  m....Mutual visibility calculation\n");
  1057.    printf ("  M....Real Time Mutual visibility.  Use the");
  1058.    printf (" following to add observers\n");
  1059.    printf ("         g-add based on grid locator\n");
  1060.    printf ("         l-add based on lat,lon,alt\n");
  1061.    printf ("         o-add based on names in qth.dat\n");
  1062.    printf ("  p....Prediction for a single satellite\n");
  1063.    printf ("  r....Real-time table of currently visible satellites,");
  1064.    printf (" additional parameter\n");
  1065.    printf ("           after r is update rate in seconds");
  1066.    printf (" (default is 1 second)\n");
  1067.    printf ("  R....Real-time table of all satellites, additional parameter");
  1068.    printf ("  after R is\n");
  1069.    printf ("           update rate in seconds (default is 1 second)\n");
  1070.    printf ("  s....Same as option R but uses a selected set of satellites\n");
  1071.    printf ("  S....Selects satellites for option s\n");
  1072.    printf ("\n");
  1073.  
  1074.    printf (" Press ENTER to continue\n");
  1075.    gets (line1);
  1076.    ClearTheScreen();
  1077.  
  1078.    printf (" EXAMPLE:\n\n");
  1079.    printf (" Real-Time Table (currently visible)...stp r 5");
  1080.    printf (" (update every 5 seconds)\n");
  1081.    printf ("\n");
  1082.    printf (" Real-Time Table (all satellites)......stp R 10");
  1083.    printf (" (update every 10 seconds)\n");
  1084.    printf ("\n");
  1085.    printf (" AOS/LOS Table (relative time).........stp a\n");
  1086.    printf ("\n");
  1087.    printf (" This version can handle up to %d Satellites",MaxSat);
  1088.    printf (" and up to %d QTH locations\n",MaxQth);
  1089.    printf ("\n");
  1090.  
  1091.    if (PlusMinus){
  1092.       printf (" Azimuth is NEGATIVE West of North and POSITIVE East of North\n");
  1093.       printf (" e.g. -90 is due west, +90 is due east\n");
  1094.    }else{
  1095.       printf (" Azimuth of 0 is north, 90 east, 180 south and 270 west\n");
  1096.    }
  1097.    printf ("\n");
  1098.  
  1099.    printf (" FindStep        = %8.2f Minutes\n",
  1100.       FindStep*1440.0);
  1101.    printf ("\n");
  1102.    printf (" FindAngleError  = %8.2f Degree(s)\n",
  1103.       FindAngleError*D2R);
  1104.    printf ("\n");
  1105.    printf (" FindStepMinimum = %8.2f Minutes\n",
  1106.       FindStepMinimum*1440.0);
  1107.    printf ("\n");
  1108.  
  1109.    printf (" PlotFlag is");
  1110.    if (PlotFlag){
  1111.       printf (" on\n");
  1112.  
  1113.       printf (" GroundTrace is");
  1114.       if (GroundTrace){
  1115.          printf (" on\n");
  1116.       }else{
  1117.          printf (" off\n");
  1118.       }
  1119.    }else{
  1120.       printf (" off\n");
  1121.    }
  1122.  
  1123.    exit(0);
  1124. }
  1125. float DayNumber(float Yinput,float Minput,float Dinput)
  1126. /* convert date to day number */
  1127. {
  1128.    float Y,M,D,DayNumb;
  1129.    float A,B;
  1130.  
  1131. /*  protect input values */
  1132.  
  1133.    Y = Yinput;
  1134.    M = Minput;
  1135.    D = Dinput;
  1136.  
  1137. /*  if january or february */
  1138.  
  1139.    if (M <= 2.0) {
  1140.       Y = Y - ONE;
  1141.       M = M + 12.0;
  1142.    }
  1143.  
  1144.    DayNumb =AINT(Y*YM) + AINT((M+ONE)*30.6001) + D;
  1145.  
  1146. /*  complete julian date calculation */
  1147.  
  1148.    A = AINT(Y/100.0);
  1149.    B = 2.0 - A + AINT(A/4.0);
  1150.    JulianDate = DayNumb + B + 1720994.5;
  1151.    
  1152.    return DayNumb - 428.0;
  1153. }
  1154. void DashLine(FILE *WhereToPrint)
  1155. /*  print dashed line */
  1156. {
  1157.    int i;
  1158.  
  1159.    for (i=0;i<=LengthOfLine;i++)fprintf (WhereToPrint,"-");
  1160.    fprintf (WhereToPrint,"\n");
  1161. }
  1162. void DataPrint(float CurrentDay)
  1163. /* print predicted data */
  1164. {
  1165.    float Altitude;
  1166.  
  1167.    PrintHeader(0,CurrentDay);
  1168.  
  1169.    fprintf (OutputFile,"%s",BLANKS);
  1170.  
  1171.    if (Step < ONE){
  1172.       Seconds = 60.0*(Minute-AINT(Minute));
  1173.       Minutes = AINT(Minute);
  1174.       fprintf (OutputFile,"%02.0f%02.0f%02.0f",
  1175.          Hour,
  1176.          Minutes,
  1177.          Seconds);
  1178.    }else{
  1179.       fprintf (OutputFile,"%02.0f%02.0f",
  1180.          Hour,
  1181.          Minute);
  1182.    }
  1183.  
  1184. /*  get mode info */
  1185.    WhatMode();
  1186.  
  1187. /*  if Omnidirectional antennas on ao-10/ao-13... then adjust squint angle */
  1188.  
  1189.    if (strncmpi(CurrentMode,"Bo",2) == 0) {
  1190.       SquintAngle = ABS(90.0 - SquintAngle);
  1191.    }
  1192.  
  1193.    if (PlotFlag){
  1194.       if (GroundTrace){
  1195.          Altitude = ZERO;
  1196.       }else{
  1197.          Altitude = AltitudeOfSatellite - RE;
  1198.       }
  1199.  
  1200.       fprintf(PlotFile1," %10.4f\t%02.0f%02.0f",
  1201.          DayNow+(Hour+Minute/60.0)/24.0-DayStart,
  1202.          Hour,
  1203.          Minute);
  1204.  
  1205.       fprintf(PlotFile1,"\t%8.1f\t%8.1f",
  1206.          Azimuth*D2R,
  1207.          Elevation*D2R);
  1208.  
  1209.       fprintf(PlotFile1,"\t%8.4f\t%8.4f\t%14.3f",
  1210.          asin(Sz/AltitudeOfSatellite)*D2R,
  1211.          ATAN2(Sy,Sx)*D2R,
  1212.          Altitude);
  1213.    }
  1214.  
  1215.    fprintf (OutputFile," %4.0f %4.0f %7.0f %4.0f",
  1216.       AINT(Azimuth*D2R+0.50),
  1217.       AINT(Elevation*D2R+0.50),
  1218.       AINT(RangeToSatellite+0.50),
  1219.       AINT(SquintAngle+0.50));
  1220.  
  1221.    Doppler(SATS[WhichSatellite].Beacon);
  1222.  
  1223. /* dump out mode */
  1224.    fprintf (OutputFile,"   %5.1f  %s",CurrentMAd,CurrentMode);
  1225.  
  1226. /* sunlit status */
  1227.    fprintf (OutputFile," %s",ECL);
  1228.  
  1229. /*  compute minimum signal level...only for ao-10 and ao-13 */
  1230.     PercivedSignalStrength(OutputFile);
  1231.  
  1232.    fprintf (OutputFile,"\n");
  1233. }
  1234. int DeltTime(float time)
  1235. /*  determine if enough time has elapsed to dump output */
  1236. {
  1237.    float dhour,dmin,dsec,tmphr,tmpmin,tmpsec;
  1238.  
  1239.    tmphr=Hour;
  1240.    tmpmin=Minute;
  1241.    tmpsec=Second;
  1242.  
  1243.    dsec  = tmpsec - SecondOld;
  1244.    if (dsec < ZERO){
  1245.       dsec += 60.0;
  1246.       tmpmin -= 1.0;
  1247.    }
  1248.  
  1249.    dmin  = tmpmin - MinuteOld;
  1250.    if (dmin < ZERO){
  1251.       dmin += 60.0;
  1252.       tmphr -= 1.0;
  1253.    }
  1254.  
  1255.    dhour = tmphr - HourOld;
  1256.    if (dhour < ZERO){
  1257.       dhour = ZERO;
  1258.    }
  1259.  
  1260.    if (dhour*3600.0+dmin*60.0+dsec >= time){
  1261.       return 1;
  1262.    }else{
  1263.       return 0;
  1264.    }
  1265. }
  1266. void DetermineAosLos(int Relative)
  1267. /*  control program to determine aos/los times */
  1268. {
  1269.    float Time,TimeStart,TimeAOS,TimeLOS,DtAOS,DtLOS;
  1270.  
  1271. /*  initailize Qth parameters */
  1272.    InitQth();
  1273.  
  1274. /*  set observer to local Qth */
  1275.  
  1276.    WhichObserver = 0;
  1277. /*  print header */
  1278.  
  1279.    if(Relative){
  1280.       printf ("        Sat       AOS            LOS\n");
  1281.    }else{
  1282.       printf ("        Sat          Day    AOS");
  1283.       printf ("               Day     LOS\n");
  1284.    }
  1285.  
  1286. /*  
  1287.   get current time....
  1288.   AOS relative to current time...
  1289.   will be zero if satellite currently visible
  1290. */
  1291.  
  1292.    UTCLOCTime();
  1293.  
  1294.    Year   = (float)GMTyear + BasicYear;
  1295.    Month  = (float)GMTmon + ONE;
  1296.    Day    = (float)GMTdaym;
  1297.  
  1298.    Hour   = (float)GMThour;
  1299.    Minute = (float)GMTmin;
  1300.    Second = (float)GMTsec;
  1301.  
  1302. /*  get day number */
  1303.  
  1304.    DayNow = DayNumber(Year,Month,Day);
  1305.  
  1306. /*  current time in hours */
  1307.  
  1308.    TimeStart = (Hour+Minute/60.0+Second/3600.0)/24.0;
  1309.  
  1310. /*  loop over each satellite and determine AOS/LOS */
  1311.  
  1312.    for (WhichSatellite=0;WhichSatellite<=NumberOfSats;WhichSatellite++){
  1313.  
  1314.       TimeNow = TimeStart;
  1315.  
  1316. /*  compute GHAE for each satellite epoch */
  1317.       InitSun();
  1318.  
  1319. /*  get current satellite location */
  1320.       SatellitePosition();
  1321.  
  1322. /*  determine azimuth and elevation from observers location */
  1323.       ObserverSatelliteGeometry();
  1324.  
  1325. /*  
  1326.    if satellite is currently up (Elevation > user value)
  1327.    then AOS problem solved (e.g. 0.0)...find LOS.
  1328. */
  1329.       if (Elevation > LOSAngle(Azimuth)){
  1330.          TimeAOS=24.0*TimeNow;
  1331.          DtAOS=ZERO;
  1332.       }else{
  1333.          TimeAOS=FindAOS();
  1334.  
  1335.          if (Relative){
  1336.             DtAOS=TimeAOS-24.0*TimeStart;
  1337.          }else{
  1338.             DtAOS=TimeAOS;
  1339.          }
  1340.       }
  1341.  
  1342. /*  find LOS */
  1343.  
  1344.       TimeLOS=FindLOS();
  1345.  
  1346.       if (Relative){
  1347.          DtLOS=TimeLOS-TimeAOS;
  1348.       }else{
  1349.          DtLOS=TimeLOS;
  1350.       }
  1351.  
  1352. /*  if relative time then dont print out day of the week */
  1353.       if (Relative)GMTdayw = -1;
  1354.  
  1355.       printf (" %12s",SATS[WhichSatellite].SatsName);
  1356.  
  1357.       if (TimeAOS < ZERO || TimeLOS < ZERO){
  1358.          printf (" Rise/Set Time > 2 Days from now!!\n");
  1359.       }else{
  1360.          printf (" ");
  1361.          HourMinuteSecond(DtAOS,stdout,GMTdayw);
  1362.  
  1363.          printf ("    ");
  1364.          HourMinuteSecond(DtLOS,stdout,GMTdayw);
  1365.          printf ("\n");
  1366.       }
  1367.    }
  1368. }
  1369. void DoCheckSum(char FileName[],int PrintFlag)
  1370. /*  driver to perform check sum computation */
  1371. {
  1372.    FILE *InFile,*OutFile;
  1373.  
  1374.    int i,nlen,NewValue,OldValue,line[80];
  1375.  
  1376.    nlen=69;
  1377.  
  1378. /*  open the input file */
  1379.  
  1380.    if ((InFile = fopen(FileName,"r")) == 0){
  1381.       printf (" Can't open file %s\n",FileName);
  1382.       exit(-1);
  1383.    }
  1384.  
  1385. /*  open the output file */
  1386.  
  1387.    if ((OutFile = fopen ("newnasa.dat","w")) == 0){
  1388.       printf (" Can't open file newnasa.dat\n");
  1389.       exit(-1);
  1390.    }
  1391.  
  1392. /*  read the file */
  1393.  
  1394.    while ((fgets(SatName,40,InFile)) != 0){
  1395.       fprintf (OutFile,"%s",SatName);
  1396.  
  1397.       NewValue = CheckSum(InFile,line);
  1398.  
  1399.       if (PrintFlag){
  1400.          OldValue = line[68] - 48;
  1401.          if (NewValue != OldValue){
  1402.             printf (" CheckSum Error for Satellite %s",SatName);
  1403.             printf ("   line 1 Old = %d  New = %d\n",
  1404.                OldValue,
  1405.                NewValue);
  1406.          }
  1407.       }
  1408.  
  1409.       line[68] = NewValue + 48;
  1410.  
  1411.       for (i=0;i<=nlen;i++){
  1412.          fprintf (OutFile,"%c",line[i]);
  1413.       }
  1414.  
  1415.       NewValue = CheckSum(InFile,line);
  1416.  
  1417.       if (PrintFlag){
  1418.          OldValue = line[68] - 48;
  1419.          if (NewValue != OldValue){
  1420.             printf (" CheckSum Error for Satellite %s",SatName);
  1421.             printf ("   line 2 Old = %d  New = %d\n",
  1422.                OldValue,
  1423.                NewValue);
  1424.          }
  1425.       }
  1426.  
  1427.       line[68] = NewValue + 48;
  1428.  
  1429.       for (i=0;i<=nlen;i++){
  1430.          fprintf (OutFile,"%c",line[i]);
  1431.       }
  1432.    }
  1433.  
  1434.    fclose (InFile);
  1435.    fclose (OutFile);
  1436. }
  1437. void Doppler(float Frequency)
  1438. /*  compute doppler shift */
  1439. {
  1440.    float doppler;
  1441.  
  1442.    doppler = -Frequency*RangeRate/299792.0;
  1443.    fprintf (OutputFile," %+8.0f",doppler);
  1444. }
  1445. void EclipseSatellitePredict(int argv,char *argc[])
  1446. /*  
  1447.   controlling routine to predict satellite eclipses
  1448.  
  1449.   really only makes sense for ao-13 and ao-10
  1450. */
  1451. {
  1452.    float CurrentDay,DayEnd,NumberOfDays;
  1453.    float ElevationLimit;
  1454.    int i,nc,RecordsWritten;
  1455.  
  1456. /*  see if user wants local eclipse's only */
  1457.  
  1458.    printf (" Do you want to determine eclipses for entire");
  1459.    printf (" orbit (press ENTER)\n");
  1460.    printf (" or determine eclipses that occur while in");
  1461.    printf (" view locally (enter a 1) :");
  1462.  
  1463.    if (GetIntNumberFromUser(&nc) == 0){
  1464.       ElevationLimit = -200.0;
  1465.    }else{
  1466.       ElevationLimit = ZERO;
  1467.    }
  1468.  
  1469.    RecordsWritten = 0;
  1470.  
  1471. /*  set observer to home Qth */
  1472.  
  1473.    WhichObserver=0;
  1474.  
  1475. /*  
  1476.   if user did not include satellite name on command line
  1477.   then show what we got and ask...
  1478. */
  1479.  
  1480.    if (argv <= 2){
  1481.       for (i=0;i<=NumberOfSats;i++){
  1482.          printf (" %s\n",SATS[i].SatsName);
  1483.       }
  1484.  
  1485.       printf (" Which Satellite from the list ?  ");
  1486.  
  1487.       gets(SatName);
  1488.    } else {
  1489.       strcpy(SatName,argc[2]);
  1490.    }
  1491.  
  1492.    nc=strlen(SatName);
  1493.  
  1494. /*  can we match it */
  1495.  
  1496.    WhichSatellite=-1;
  1497.  
  1498.    for (i=0;i<=NumberOfSats;i++){
  1499.       if (strnicmp(SatName,SATS[i].SatsName,nc) == 0){
  1500.          if (WhichSatellite == -1){
  1501.             WhichSatellite = i;
  1502.          }else{
  1503.             printf (" Multiple Names.  Please be more specific \n");
  1504.             exit(-1);
  1505.          }
  1506.       }
  1507.    }
  1508.  
  1509.    if (WhichSatellite < 0){
  1510.       printf (" Could Not Match Satellite %s in List\n",SatName);
  1511.       exit (-1);
  1512.    }else{
  1513.       printf (" Matched your input as %s\n",SATS[WhichSatellite].SatsName);
  1514.    }
  1515.  
  1516. /*  
  1517.   get local date to check input (currently only warns user of 
  1518.   back propagation) and possibly for use as input
  1519. */
  1520.    UTCLOCTime();
  1521.  
  1522.    Month = (float)GMTmon + ONE;
  1523.    Day   = (float)GMTdaym;
  1524.    Year  = (float)GMTyear + BasicYear;
  1525.  
  1526. /*  get current day number */
  1527.    CurrentDay=DayNumber(Year,Month,Day);
  1528.  
  1529.    printf (" Enter Start Day (month day year) (<cr> for current date) : ");
  1530.  
  1531. /*  allow user option to hit return to get current date */
  1532.    strcpy(line1,GetStringFromUser(&nc));
  1533.  
  1534. /*  if nc > 0 then user input...otherwise use current information */
  1535.  
  1536.    if (nc > 0){
  1537.       sscanf(line1,"%f%f%f",&Month,&Day,&Year);
  1538.    }else{
  1539.       Month = (float)GMTmon + ONE;
  1540.       Day   = (float)GMTdaym;
  1541.       Year  = (float)GMTyear + BasicYear;
  1542.    }
  1543.  
  1544. /*  check for back propagation times */
  1545.  
  1546.    if (Month < (float)LOCmon){
  1547.       printf (" ********* Input Month < Current Month\n\n");
  1548.    }
  1549.    if (Day < (float)LOCdaym && Month <= (float)LOCmon){
  1550.       printf (" ********* Input Day < Current Day\n\n");
  1551.    }
  1552.    if (Year < BasicYear +(float)LOCyear){
  1553.       printf (" ********* Input Year < Current Year\n");
  1554.       printf ("           MA calculations will be wrong\n\n");
  1555.    }
  1556.  
  1557.    printf (" Number of Days for Prediction: ");
  1558.    NumberOfDays=GetFloatNumberFromUser(&nc);
  1559.  
  1560.    printf (" Step Size (minutes....input negative value for seconds): ");
  1561.    Step=GetFloatNumberFromUser(&nc);
  1562.  
  1563. /* 
  1564.    if Step is negative then user has input time step in seconds...
  1565.    convert to minutes 
  1566. */
  1567.  
  1568.    if (Step < ZERO){
  1569.       Step = ABS(Step)/60.0;
  1570.    }
  1571.  
  1572.    printf (" FileName to write output to (cr for screen): ");
  1573.  
  1574. /* read filename here */
  1575.  
  1576.    strcpy(FileName,GetStringFromUser(&nc));
  1577.  
  1578.    if (nc > 0){
  1579.       if ((OutputFile = fopen(FileName,"w")) == 0){
  1580.          printf ("Can't open file %s\n",FileName);
  1581.          exit(-1);
  1582.       }
  1583.    }else{
  1584.       OutputFile = stdout;
  1585.    }
  1586.  
  1587.    DayStart = DayNumber(Year,Month,Day);
  1588.    DayEnd   = DayStart + NumberOfDays - ONE;
  1589.  
  1590. /*  initailize Qth parameters */
  1591.  
  1592.    InitQth();
  1593.  
  1594. /*  initailize sun parameters */
  1595.  
  1596.    InitSun();
  1597.  
  1598. /*  initailize antenna parameters */
  1599.  
  1600.    InitAntenna();
  1601.  
  1602. /*  loop over time */
  1603.  
  1604.    RunningCounter    = 0;
  1605.    RunningCounterOld = 0;
  1606.  
  1607.    for (DayNow=DayStart;DayNow<=DayEnd;DayNow++){
  1608.       for (Hour=ZERO;Hour<=23.0;Hour++){
  1609.          for (Minute=ZERO;Minute<=60.0-Step;Minute=Minute+Step){
  1610.  
  1611.             RunningCounter++;
  1612.  
  1613.             TimeNow = (Hour+Minute/60.0)/24.0;
  1614.  
  1615. /*  get satellite location */
  1616.             SatellitePosition();
  1617.  
  1618. /*  determine azimuth and elevation from observers location */
  1619.             ObserverSatelliteGeometry();
  1620.  
  1621. /*  
  1622.    since we want total eclipse information check constanly
  1623.  
  1624.    compute sun illumination 
  1625.    squint angle (antenna pointing)
  1626. */
  1627.             SunData();
  1628.  
  1629. /*  if eclipsed then dump time */
  1630.  
  1631.             if ((strcmp(ECL," ECL ") == 0) && Elevation > ElevationLimit){
  1632.                SquintRangeRate();
  1633.                DataPrint(CurrentDay);
  1634.                RecordsWritten++;
  1635.             }
  1636.          }
  1637.       }
  1638.    }
  1639.  
  1640.    if (RecordsWritten == 0){
  1641.       printf (" No eclipses found during the time interval!!\n");
  1642.    }
  1643.  
  1644.    fclose (OutputFile);
  1645. }
  1646. void EditKepler(void)
  1647. /*  allow user to edit satellite parameters in keplerian format*/
  1648. {
  1649.    FILE *EditFile;
  1650.  
  1651.    char dummy[2];
  1652.    int i,ndigit,ThisSat;
  1653.    float EpochDay,tempvar;
  1654.  
  1655.    for (i=0;i<=NumberOfSats;i++){
  1656.       printf (" %4d %s\n",i,SATS[i].SatsName);
  1657.    }
  1658.  
  1659.    printf (" Which sat (0-%d) ",NumberOfSats);
  1660.    ThisSat=InputCheck(0,NumberOfSats);
  1661.  
  1662.    if (ThisSat <= NumberOfSats){
  1663.       EpochDay=SATS[ThisSat].DE-
  1664.               DayNumber(SATS[ThisSat].YE,ONE,ZERO)+
  1665.               SATS[ThisSat].TE;
  1666.  
  1667.       printf (" Satellite      : %s \n",SATS[ThisSat].SatsName);
  1668.       printf (" Catalog Number : %15d \n",SATS[ThisSat].SatNum);
  1669.       printf (" Element Set    : %15d \n",SATS[ThisSat].EleSet);
  1670.       printf (" Epoch (year)   : %15.10f \n",SATS[ThisSat].YE-BasicYear);
  1671.       printf (" Epoch (day)    : %15.10f \n",EpochDay);
  1672.       printf (" Inclination    : %15.10f \n",SATS[ThisSat].IN*D2R);
  1673.       printf (" R.A.A.N        : %15.10f \n",SATS[ThisSat].RA*D2R);
  1674.       printf (" Eccentricity   : %15.10f \n",SATS[ThisSat].EC);
  1675.       printf (" Arg of Perigee : %15.10f \n",SATS[ThisSat].WP*D2R);
  1676.       printf (" Mean Anomaly   : %15.10f \n",SATS[ThisSat].MA*D2R);
  1677.       printf (" Mean Motion    : %15.10f \n",SATS[ThisSat].MM/TWOPI);
  1678.       printf (" Decay          : %15.10f \n",SATS[ThisSat].M2/TWOPI);
  1679.       printf (" Orbit #        : %15.10f \n",SATS[ThisSat].RV);
  1680.  
  1681.       printf ("\n");
  1682.       printf (" press return to leave value alone..");
  1683.       printf ("  enter new value to change\n\n");
  1684.  
  1685. /*  get carriage return from input line */
  1686.  
  1687.       gets(dummy);
  1688.  
  1689.       printf (" Element Set    : %15d ",SATS[ThisSat].EleSet);
  1690.       tempvar=GetFloatNumberFromUser(&ndigit);
  1691.       if (ndigit > 0){
  1692.          SATS[ThisSat].EleSet=(int)tempvar;
  1693.          printf (" new value %15.10f\n",tempvar);
  1694.       }else{
  1695.          SATS[ThisSat].EleSet=SATS[ThisSat].EleSet;
  1696.       }
  1697.  
  1698.       printf (" Epoch (year)   : %15.10f ",SATS[ThisSat].YE-BasicYear);
  1699.       tempvar=GetFloatNumberFromUser(&ndigit);
  1700.       if (ndigit > 0){
  1701.          SATS[ThisSat].YE=tempvar+BasicYear;
  1702.          printf (" new value %15.10f\n",tempvar);
  1703.       }else{
  1704.          SATS[ThisSat].YE=SATS[ThisSat].YE;
  1705.       }
  1706.  
  1707.       printf (" Epoch (day)    : %15.10f ",EpochDay);
  1708.       tempvar=GetFloatNumberFromUser(&ndigit);
  1709.       if (ndigit > 0){
  1710.          SATS[ThisSat].TE=tempvar;
  1711.          printf (" new value %15.10f\n",tempvar);
  1712.       }else{
  1713.          SATS[ThisSat].TE=EpochDay;
  1714.       }
  1715.  
  1716.       printf (" Inclination    : %15.10f ",SATS[ThisSat].IN*D2R);
  1717.       tempvar=GetFloatNumberFromUser(&ndigit);
  1718.       if (ndigit > 0){
  1719.          SATS[ThisSat].IN=tempvar;
  1720.          printf (" new value %15.10f\n",tempvar);
  1721.       }else{
  1722.          SATS[ThisSat].IN=SATS[ThisSat].IN;
  1723.       }
  1724.  
  1725.       printf (" R.A.A.N        : %15.10f ",SATS[ThisSat].RA*D2R);
  1726.       tempvar=GetFloatNumberFromUser(&ndigit);
  1727.       if (ndigit > 0){
  1728.          SATS[ThisSat].RA=tempvar;
  1729.          printf (" new value %15.10f\n",tempvar);
  1730.       }else{
  1731.          SATS[ThisSat].RA=SATS[ThisSat].RA;
  1732.       }
  1733.  
  1734.       printf (" Eccentricity   : %15.10f ",SATS[ThisSat].EC);
  1735.       tempvar=GetFloatNumberFromUser(&ndigit);
  1736.       if (ndigit > 0){
  1737.          SATS[ThisSat].EC=tempvar;
  1738.          printf (" new value %15.10f\n",tempvar);
  1739.       }else{
  1740.          SATS[ThisSat].EC=SATS[ThisSat].EC;
  1741.       }
  1742.  
  1743.       printf (" Arg of Perigee : %15.10f ",SATS[ThisSat].WP*D2R);
  1744.       tempvar=GetFloatNumberFromUser(&ndigit);
  1745.       if (ndigit > 0){
  1746.          SATS[ThisSat].WP=tempvar;
  1747.          printf (" new value %15.10f\n",tempvar);
  1748.       }else{
  1749.          SATS[ThisSat].WP=SATS[ThisSat].WP;
  1750.       }
  1751.  
  1752.       printf (" Mean Anomaly   : %15.10f ",SATS[ThisSat].MA*D2R);
  1753.       tempvar=GetFloatNumberFromUser(&ndigit);
  1754.       if (ndigit > 0){
  1755.          SATS[ThisSat].MA=tempvar;
  1756.          printf (" new value %15.10f\n",tempvar);
  1757.       }else{
  1758.          SATS[ThisSat].MA=SATS[ThisSat].MA;
  1759.       }
  1760.  
  1761.       printf (" Mean Motion    : %15.10f ",SATS[ThisSat].MM/TWOPI);
  1762.       tempvar=GetFloatNumberFromUser(&ndigit);
  1763.       if (ndigit > 0){
  1764.          SATS[ThisSat].MM=tempvar;
  1765.          printf (" new value %15.10f\n",tempvar);
  1766.       }else{
  1767.          SATS[ThisSat].MM=SATS[ThisSat].MM;
  1768.       }
  1769.  
  1770.       printf (" Decay          : %15.10f ",SATS[ThisSat].M2/TWOPI);
  1771.       tempvar=GetFloatNumberFromUser(&ndigit);
  1772.       if (ndigit > 0){
  1773.          SATS[ThisSat].M2=tempvar;
  1774.          printf (" new value %15.10f\n",tempvar);
  1775.       }else{
  1776.          SATS[ThisSat].M2=SATS[ThisSat].M2;
  1777.       }
  1778.  
  1779.       printf (" Orbit #        : %15.10f ",SATS[ThisSat].RV);
  1780.       tempvar=GetFloatNumberFromUser(&ndigit);
  1781.       if (ndigit > 0){
  1782.          SATS[ThisSat].RV=tempvar;
  1783.          printf (" new value %15.10f\n",tempvar);
  1784.       }else{
  1785.          SATS[ThisSat].RV=SATS[ThisSat].RV;
  1786.       }
  1787.    }
  1788.  
  1789. /*  if here then open new file and dump modified data......kepler format */
  1790.  
  1791.    if((EditFile=fopen("kepler.tmp","w")) == 0){
  1792.       printf (" cant open kepler.tmp\n");
  1793.       exit(-1);
  1794.    }
  1795.  
  1796.    for (i=0;i<=NumberOfSats;i++){
  1797.       fprintf (EditFile,"Satellite: %s\n",SATS[i].SatsName);
  1798.       fprintf (EditFile,"Catalog number: %d\n",SATS[i].SatNum);
  1799.  
  1800.       tempvar=SATS[i].YE-BasicYear;
  1801.  
  1802. /*  
  1803.    if satellite we modified then dump new info...
  1804.    otherwise generate correct date
  1805. */
  1806.       if (i == ThisSat){
  1807.          EpochDay=SATS[i].TE;
  1808.       }else{
  1809.          EpochDay=SATS[i].DE-DayNumber(SATS[i].YE,ONE,ZERO)+SATS[i].TE;
  1810.       }
  1811.  
  1812.       fprintf (EditFile,"Epoch time:     %02.0f%9.6f\n",tempvar,EpochDay);
  1813.       fprintf (EditFile,"Element set:    %d\n",SATS[i].EleSet);
  1814.       fprintf (EditFile,"Inclination:    %12.6f\n",SATS[i].IN*D2R);
  1815.       fprintf (EditFile,"RA of node:     %12.6f\n",SATS[i].RA*D2R);
  1816.       fprintf (EditFile,"Eccentricity:   %12.6f\n",SATS[i].EC);
  1817.       fprintf (EditFile,"Arg of perigee: %12.6f\n",SATS[i].WP*D2R);
  1818.       fprintf (EditFile,"Mean anomaly:   %12.6f\n",SATS[i].MA*D2R);
  1819.       fprintf (EditFile,"Mean motion:    %12.6f\n",SATS[i].MM/TWOPI);
  1820.       fprintf (EditFile,"Decay rate:     %12.10f\n",SATS[i].M2/TWOPI);
  1821.       fprintf (EditFile,"Epoch rev:      %7.0f\n",SATS[i].RV);
  1822.       fprintf (EditFile,"\n");
  1823.    }
  1824. }
  1825. void EditNasa(void)
  1826. /*  edit data in nasa 2 line format */
  1827. {
  1828.    FILE *EditFile;
  1829.  
  1830.    char dummy[2];
  1831.    int ndigit,ThisSat,i,nc,idum,YearLau,LauNum;
  1832.    unsigned long int EpochRev;
  1833.    float EpochYear,EpochDay,DecayRate,MeanMotion,MeanAnomaly;
  1834.    float Inclination,RAofAN,Eccentricity,ArgPerigee,tempvar;
  1835.  
  1836.    printf ("\a\a\a\a\n");
  1837.    printf (" WARNING..... STP does not keep all the data stored in\n");
  1838.    printf ("              a 2-line element file.  It only keeps what\n");
  1839.    printf ("              is necessary for this formulation......\n");
  1840.    printf ("\n");
  1841.  
  1842.    for (i=0;i<=NumberOfSats;i++){
  1843.       printf (" %4d %s\n",i,SATS[i].SatsName);
  1844.    }
  1845.  
  1846.    printf (" Which sat (0-%d) ",NumberOfSats);
  1847.    ThisSat=InputCheck(0,NumberOfSats);
  1848.  
  1849.    if ((EditFile = fopen("nasa.tmp","w")) == 0){
  1850.       printf (" cant open nasa.tmp\n");
  1851.       exit(-1);
  1852.    }
  1853.  
  1854.    EpochDay    = SATS[ThisSat].DE 
  1855.                 -DayNumber(SATS[ThisSat].YE,ONE,ZERO)
  1856.                 +SATS[ThisSat].TE;
  1857.  
  1858.    printf (" Satellite      : %s \n",SATS[ThisSat].SatsName);
  1859.    printf (" Catalog Number : %15d \n",SATS[ThisSat].SatNum);
  1860.    printf (" Element Set    : %15d \n",SATS[ThisSat].EleSet);
  1861.    printf (" Epoch (year)   : %15.10f \n",SATS[ThisSat].YE-BasicYear);
  1862.    printf (" Epoch (day)    : %15.10f \n",EpochDay);
  1863.    printf (" Inclination    : %15.10f \n",SATS[ThisSat].IN*D2R);
  1864.    printf (" R.A.A.N        : %15.10f \n",SATS[ThisSat].RA*D2R);
  1865.    printf (" Eccentricity   : %15.10f \n",SATS[ThisSat].EC);
  1866.    printf (" Arg of Perigee : %15.10f \n",SATS[ThisSat].WP*D2R);
  1867.    printf (" Mean Anomaly   : %15.10f \n",SATS[ThisSat].MA*D2R);
  1868.    printf (" Mean Motion    : %15.10f \n",SATS[ThisSat].MM/TWOPI);
  1869.    printf (" Decay          : %15.10f \n",SATS[ThisSat].M2/TWOPI);
  1870.    printf (" Orbit #        : %15.10f \n",SATS[ThisSat].RV);
  1871.  
  1872.    printf ("\n");
  1873.    printf (" press return to leave value alone..");  
  1874.    printf ("  enter new value to change\n\n");
  1875.  
  1876. /*  get carriage return from input line */
  1877.  
  1878.    gets(dummy);
  1879.  
  1880.    printf (" Element Set    : %15d ",SATS[ThisSat].EleSet);
  1881.    tempvar=GetFloatNumberFromUser(&ndigit);
  1882.    if (ndigit > 0){
  1883.       SATS[ThisSat].EleSet=(int)tempvar;
  1884.       printf (" new value %15.10f\n",tempvar);
  1885.    }else{
  1886.       SATS[ThisSat].EleSet=SATS[ThisSat].EleSet;
  1887.    }
  1888.  
  1889.    printf (" Epoch (year)   : %15.10f ",SATS[ThisSat].YE-BasicYear);
  1890.    tempvar=GetFloatNumberFromUser(&ndigit);
  1891.    if (ndigit > 0){
  1892.       SATS[ThisSat].YE=tempvar+BasicYear;
  1893.       printf (" new value %15.10f\n",tempvar);
  1894.    }else{
  1895.       SATS[ThisSat].YE=SATS[ThisSat].YE;
  1896.    }
  1897.  
  1898.    printf (" Epoch (day)    : %15.10f ",EpochDay);
  1899.    tempvar=GetFloatNumberFromUser(&ndigit);
  1900.    if (ndigit > 0){
  1901.       SATS[ThisSat].TE=tempvar;
  1902.       printf (" new value %15.10f\n",tempvar);
  1903.    }else{
  1904.       SATS[ThisSat].TE=EpochDay;
  1905.    }
  1906.  
  1907.    printf (" Inclination    : %15.10f ",SATS[ThisSat].IN*D2R);
  1908.    tempvar=GetFloatNumberFromUser(&ndigit);
  1909.    if (ndigit > 0){
  1910.       SATS[ThisSat].IN=tempvar;
  1911.       printf (" new value %15.10f\n",tempvar);
  1912.    }else{
  1913.       SATS[ThisSat].IN=SATS[ThisSat].IN;
  1914.    }
  1915.  
  1916.    printf (" R.A.A.N        : %15.10f ",SATS[ThisSat].RA*D2R);
  1917.    tempvar=GetFloatNumberFromUser(&ndigit);
  1918.    if (ndigit > 0){
  1919.       SATS[ThisSat].RA=tempvar;
  1920.       printf (" new value %15.10f\n",tempvar);
  1921.    }else{
  1922.       SATS[ThisSat].RA=SATS[ThisSat].RA;
  1923.    }
  1924.  
  1925.    printf (" Eccentricity   : %15.10f ",SATS[ThisSat].EC);
  1926.    tempvar=GetFloatNumberFromUser(&ndigit);
  1927.    if (ndigit > 0){
  1928.       SATS[ThisSat].EC=tempvar;
  1929.       printf (" new value %15.10f\n",tempvar);
  1930.    }else{
  1931.       SATS[ThisSat].EC=SATS[ThisSat].EC;
  1932.    }
  1933.  
  1934.    printf (" Arg of Perigee : %15.10f ",SATS[ThisSat].WP*D2R);
  1935.    tempvar=GetFloatNumberFromUser(&ndigit);
  1936.    if (ndigit > 0){
  1937.       SATS[ThisSat].WP=tempvar;
  1938.       printf (" new value %15.10f\n",tempvar);
  1939.    }else{
  1940.       SATS[ThisSat].WP=SATS[ThisSat].WP;
  1941.    }
  1942.  
  1943.    printf (" Mean Anomaly   : %15.10f ",SATS[ThisSat].MA*D2R);
  1944.    tempvar=GetFloatNumberFromUser(&ndigit);
  1945.    if (ndigit > 0){
  1946.       SATS[ThisSat].MA=tempvar;
  1947.       printf (" new value %15.10f\n",tempvar);
  1948.    }else{
  1949.       SATS[ThisSat].MA=SATS[ThisSat].MA;
  1950.    }
  1951.  
  1952.    printf (" Mean Motion    : %15.10f ",SATS[ThisSat].MM/TWOPI);
  1953.    tempvar=GetFloatNumberFromUser(&ndigit);
  1954.    if (ndigit > 0){
  1955.       SATS[ThisSat].MM=tempvar;
  1956.       printf (" new value %15.10f\n",tempvar);
  1957.    }else{
  1958.       SATS[ThisSat].MM=SATS[ThisSat].MM;
  1959.    }
  1960.  
  1961.    printf (" Decay          : %15.10f ",SATS[ThisSat].M2/TWOPI);
  1962.    tempvar=GetFloatNumberFromUser(&ndigit);
  1963.    if (ndigit > 0){
  1964.       SATS[ThisSat].M2=tempvar;
  1965.       printf (" new value %15.10f\n",tempvar);
  1966.    }else{
  1967.       SATS[ThisSat].M2=SATS[ThisSat].M2;
  1968.    }
  1969.  
  1970.    printf (" Orbit #        : %15.10f ",SATS[ThisSat].RV);
  1971.    tempvar=GetFloatNumberFromUser(&ndigit);
  1972.    if (ndigit > 0){
  1973.       SATS[ThisSat].RV=tempvar;
  1974.       printf (" new value %15.10f\n",tempvar);
  1975.    }else{
  1976.       SATS[ThisSat].RV=SATS[ThisSat].RV;
  1977.    }
  1978.  
  1979. /*  dump out new data */
  1980.  
  1981.    for (i=0;i<=NumberOfSats;i++){
  1982.  
  1983. /*
  1984.    if satellite we modified then dump new info...
  1985.    otherwise generate correct date
  1986. */
  1987.       if (i == ThisSat){
  1988.          EpochDay = SATS[i].TE;
  1989.       }else{
  1990.          EpochDay = SATS[i].DE
  1991.                    -DayNumber(SATS[i].YE,ONE,ZERO)
  1992.                    +SATS[i].TE;
  1993.       }     
  1994.  
  1995.       EpochRev    = (int)SATS[i].RV;
  1996.       EpochYear   = SATS[i].YE - BasicYear;
  1997.       DecayRate   = SATS[i].M2/TWOPI;
  1998.       MeanMotion  = SATS[i].MM/TWOPI;
  1999.       MeanAnomaly = SATS[i].MA*D2R;
  2000.       Inclination = SATS[i].IN*D2R;
  2001.       RAofAN      = SATS[i].RA*D2R;
  2002.       Eccentricity= SATS[i].EC;
  2003.       ArgPerigee  = SATS[i].WP*D2R;
  2004.  
  2005.  
  2006. /*  satellite name */
  2007.  
  2008.       fprintf (EditFile,"%s\n",SATS[i].SatsName);
  2009.  
  2010. /*  1st line */
  2011.  
  2012.       fprintf (EditFile,"1 ");
  2013.       fprintf (EditFile,"%5dU",
  2014.          SATS[i].SatNum);
  2015.       fprintf (EditFile," 00000A");
  2016.       fprintf (EditFile,"   %2.0f%11.8f",
  2017.          EpochYear,
  2018.          EpochDay);
  2019.       if (DecayRate < ZERO){
  2020.          fprintf (EditFile," -.%08.0f",
  2021.             ABS(DecayRate)*1.0e8);
  2022.       }else{
  2023.          fprintf (EditFile,"  .%08.0f",
  2024.             DecayRate*1.0e8);
  2025.       }
  2026.       fprintf (EditFile,"  00000-0");
  2027.       fprintf (EditFile,"  10000-3 0");
  2028.       fprintf (EditFile,"% 5d",
  2029.          SATS[i].EleSet);
  2030.       fprintf (EditFile,"0");
  2031.       fprintf (EditFile,"\n");
  2032.  
  2033. /*  2nd line */
  2034.  
  2035.       fprintf (EditFile,"2 ");
  2036.       fprintf (EditFile,"%5d",
  2037.          SATS[i].SatNum);
  2038.       fprintf (EditFile," %8.4f",
  2039.          Inclination);
  2040.       fprintf (EditFile," %8.4f",
  2041.          RAofAN);
  2042.       fprintf (EditFile," %07.0f",
  2043.          Eccentricity*1.0e7);
  2044.       fprintf (EditFile," %8.4f",
  2045.          ArgPerigee);
  2046.       fprintf (EditFile," %8.4f",
  2047.          MeanAnomaly);
  2048.       fprintf (EditFile,"% 12.8f",
  2049.          MeanMotion);
  2050.       fprintf (EditFile,"%5d",
  2051.          EpochRev*10);
  2052.       fprintf (EditFile,"0");
  2053.       fprintf (EditFile,"\n");
  2054.    }
  2055.  
  2056. /*  close file and reread and write checksums */
  2057.    fclose (EditFile);
  2058.  
  2059.    DoCheckSum("new.tmp",0);
  2060. }
  2061. float FindAOS(void)
  2062. /* 
  2063.     find acquisition time...
  2064.  
  2065.     StepSize algorithm based on work by NZ3F and
  2066.     WN2A (author of SatSked)
  2067.  
  2068.     use current elevation angle to determine stepsize
  2069.     big angle....big stepsize
  2070. */
  2071. {
  2072.    float Stepp,TotalElapsedTime;
  2073.    float xneg,yneg;
  2074.  
  2075.    Stepp = FindStep*ABS(Elevation*Elevation);
  2076.  
  2077. /* 
  2078.    step time forward and determine elevation angle.
  2079.    when we find a elevation angle > requirment then
  2080.    we linear interpolate between last negative answer
  2081.    and this one to get the answer
  2082. */
  2083.  
  2084.    TotalElapsedTime = ZERO;
  2085.  
  2086.    while (Elevation < ZERO){
  2087.       yneg=Elevation;
  2088.       xneg=TimeNow;
  2089.  
  2090.       TimeNow += Stepp;
  2091.  
  2092.       TotalElapsedTime += Stepp;
  2093.       if (TotalElapsedTime > 2.0){
  2094.          return -ONE;
  2095.       }
  2096.  
  2097.       SatellitePosition();
  2098.       ObserverSatelliteGeometry();
  2099.  
  2100.       Stepp = FindStep*ABS(Elevation*Elevation);
  2101.       if (Stepp < FindStepMinimum) Stepp = FindStepMinimum;
  2102.    }
  2103.  
  2104. /*  if within error then return */
  2105.     
  2106.     if (ABS(Elevation) < FindAngleError){
  2107.        return 24.0*TimeNow;
  2108.     }
  2109.  
  2110. /*  get better answer */
  2111.  
  2112.    TimeNow = xneg;
  2113.    Elevation = yneg;
  2114.  
  2115.    Stepp /= 10.0;
  2116.  
  2117.    while (Elevation < ZERO){
  2118.       yneg=Elevation;
  2119.  
  2120.       TimeNow += Stepp;
  2121.  
  2122.       SatellitePosition();
  2123.       ObserverSatelliteGeometry();
  2124.  
  2125.       Stepp = FindStep*ABS(Elevation*Elevation);
  2126.       if (Stepp < FindStepMinimum) Stepp = FindStepMinimum;
  2127.    }
  2128.  
  2129.    return 24.0*TimeNow;
  2130. }
  2131. float FindLOS()
  2132. /*  find out when a satellite goes below horizion limit */
  2133. {
  2134.    float Stepp,TotalElapsedTime;
  2135.   
  2136.    float xpos,ypos;
  2137.  
  2138.    Stepp = FindStep*ABS(Elevation*Elevation);
  2139. /* 
  2140.    step time forward and determine elevation angle.
  2141.    when we find a elevation angle < requirment then
  2142.    we linear interpolate between last positive answer
  2143.    and this one to get the answer
  2144. */
  2145.  
  2146.    TotalElapsedTime = ZERO;
  2147.  
  2148.    while (Elevation > ZERO){
  2149.  
  2150.       xpos=TimeNow;
  2151.       ypos=Elevation;
  2152.  
  2153.       TimeNow += Stepp;
  2154.  
  2155.       TotalElapsedTime += Stepp;
  2156.       if (TotalElapsedTime > 2.0){
  2157.          return -ONE;
  2158.       }
  2159.  
  2160.       SatellitePosition();
  2161.       ObserverSatelliteGeometry();
  2162.  
  2163.       Stepp = FindStep*ABS(Elevation*Elevation);
  2164.       if (Stepp < FindStepMinimum) Stepp = FindStepMinimum;
  2165.    }
  2166.  
  2167. /*  if within limit then return */
  2168.     
  2169.     if (ABS(Elevation) < FindAngleError){
  2170.        return 24.0*TimeNow;
  2171.     }
  2172.  
  2173. /*  get better answer */
  2174.  
  2175.    TimeNow = xpos;
  2176.    Elevation = ypos;
  2177.  
  2178.    Stepp /= 10.0;
  2179.  
  2180.    while (Elevation > ZERO){
  2181.       ypos=Elevation;
  2182.  
  2183.       TimeNow += Stepp;
  2184.  
  2185.       SatellitePosition();
  2186.       ObserverSatelliteGeometry();
  2187.  
  2188.       Stepp = FindStep*ABS(Elevation*Elevation);
  2189.       if (Stepp < FindStepMinimum) Stepp = FindStepMinimum;
  2190.    }
  2191.  
  2192.    return 24.0*TimeNow;
  2193. }
  2194. float FNatn(float Y,float X)
  2195. /*  G3RUH atan function....*/
  2196. {
  2197.    float A;
  2198.  
  2199.    if (X == ZERO){
  2200.       A = PIOVR2*SGN(Y);
  2201.    }else{
  2202.       A = atan(Y/X);
  2203.    }
  2204.    if (X < ZERO) A += PI;
  2205.    if (A < ZERO) A += TWOPI;
  2206.  
  2207.    return A;
  2208. }
  2209. void FootPrintInit(void)
  2210. {
  2211.    int i;
  2212.    float A;
  2213.  
  2214.    NumberOfFootprint = 40;
  2215.  
  2216.    if (NumberOfFootprint > MaxFootprint){
  2217.       printf (" NumberofFootprint > MaxFootprint %d %d\n",
  2218.          NumberOfFootprint,
  2219.          MaxFootprint);
  2220.       exit(-1);
  2221.    }
  2222.  
  2223.    for (i=0;i<=NumberOfFootprint;i++){
  2224.       A = TWOPI*(float)i/(float)NumberOfFootprint;
  2225.       SinOfA[i] = sin(A);
  2226.       CosOfA[i] = cos(A);
  2227.    }
  2228. }
  2229. void FootPrint(void)
  2230. /*  determine lat,lon on ground that defines footprint of satellite */
  2231. {
  2232.    int i;
  2233.    float RadiusOfFootprint,SubSatLatitude,SubSatLongitude;
  2234.    float CosLat,CosLon,SinLat,SinLon,CosRad,SinRad,X,Y,Z,x,y,z;
  2235.  
  2236. /*  maximum distance from subsatellite point that can view satellite */
  2237.  
  2238.    RadiusOfFootprint = acos(RE/AltitudeOfSatellite);
  2239.  
  2240. /*  compute subsatellite point at current time */
  2241.  
  2242.    SubSatLatitude  = asin(Sz/AltitudeOfSatellite);
  2243.    SubSatLongitude = ATAN2(Sy,Sx);
  2244.  
  2245. /*  pre compute some constants */
  2246.  
  2247.    CosLat = cos(SubSatLatitude);
  2248.    CosLon = cos(SubSatLongitude);
  2249.    SinLat = sin(SubSatLatitude);
  2250.    SinLon = sin(SubSatLongitude);
  2251.  
  2252.    CosRad = cos(RadiusOfFootprint);
  2253.    SinRad = sin(RadiusOfFootprint);
  2254.  
  2255. /*  determine the points */
  2256.  
  2257.    for (i=0;i<=NumberOfFootprint;i++){
  2258.  
  2259. /*  point centered at 0,0 */
  2260.  
  2261.       X = CosRad;
  2262.       Y = SinRad*SinOfA[i];
  2263.       Z = SinRad*CosOfA[i];
  2264.  
  2265. /*  rotate up by latitude */
  2266.  
  2267.       x = X*CosLat - Z*SinLat;
  2268.       y = Y;
  2269.       z = X*SinLat + Z*CosLat;
  2270.  
  2271. /*  rotate around by longitude */
  2272.  
  2273.       X = x*CosLon - y*SinLon;
  2274.       Y = x*SinLon + y*CosLon;
  2275.       Z = z;
  2276.  
  2277. /*  convert to altitude and longitude */
  2278.  
  2279.       FootprintLat[i] = asin(Z)*D2R;
  2280.       FootprintLon[i] = ATAN2(Y,X)*D2R;
  2281.  
  2282.       fprintf (PlotFile2," %10.3f %10.3f %4.0f %02.0f%02.0f\n",
  2283.          FootprintLat[i],
  2284.          FootprintLon[i],
  2285.          ZERO,
  2286.          Hour,Minute);
  2287.    }
  2288. }
  2289. float GetFloatNumberFromUser(int *m)
  2290. /*
  2291.   return value of number from line....if none entered return 0.
  2292.   also return number of values on line (m)
  2293.  
  2294.   method suggested my eric odell
  2295. */
  2296. {
  2297.    char numbers[40];
  2298.    int ich,n;
  2299.  
  2300.    n=0;
  2301.  
  2302. /*  
  2303.   read numbers from line and store as characters
  2304.   when you get to carriage return then finished
  2305. */
  2306.  
  2307.    while ((ich=getchar()) !='\n'){
  2308.       numbers[n]=ich;
  2309.       n++;
  2310.    }
  2311.  
  2312. /*  put end of line mark for atof function */
  2313.  
  2314.    numbers[n]='\0';
  2315.  
  2316.    *m=n;
  2317.  
  2318. /*  if n > 0 then return atof(numbers), otherwise return 0.0 */
  2319.  
  2320.    return (n) ? atof(numbers): 0.0;
  2321. }
  2322. int GetIntNumberFromUser(int *m)
  2323. /*
  2324.   return int value of number from line....if none entered return 0.
  2325.   also return number of values on line (m)
  2326. */
  2327. {
  2328.    char numbers[40];
  2329.    int ich,n;
  2330.  
  2331.    n=0;
  2332.  
  2333. /*
  2334.   read numbers from line and store as characters
  2335.   when you get to carriage return then finished
  2336. */
  2337.  
  2338.    while ((ich=getchar()) !='\n'){
  2339.       numbers[n]=ich;
  2340.       n++;
  2341.    }
  2342.  
  2343. /*  put end of line mark for atoi function */
  2344.  
  2345.    numbers[n]='\0';
  2346.  
  2347.    *m=n;
  2348.  
  2349. /*  if n > 0 then return atoi(numbers), otherwise return 0 */
  2350.  
  2351.    return (n) ? atoi(numbers): 0;
  2352. }
  2353. char *GetStringFromUser(int *m)
  2354. /*
  2355.   return string from line....if none entered return 0.
  2356.   also return number of letters on line (m)
  2357. */
  2358. {
  2359.    char garbage[132],numbers[132];
  2360.    int ich,n,l;
  2361.  
  2362.    n=0;
  2363.  
  2364. /*
  2365.   read letters from line and store.
  2366.   when you get to carriage return then finished
  2367. */
  2368.  
  2369.    while ((ich=getchar()) !='\n'){
  2370.       garbage[n]=ich;
  2371.       n++;
  2372.    }
  2373.  
  2374.    garbage[n]='\0';
  2375.    strcpy(numbers,garbage);
  2376.  
  2377. /*  put end of line mark for atof function */
  2378.  
  2379.    numbers[n]='\0';                           
  2380.  
  2381.    *m=n;            
  2382.  
  2383. /*  if n > 0 then return string, otherwise return NULL */
  2384.  
  2385.    return (n) ? numbers : BLANK;                          
  2386. }      
  2387. int GridLocator(char *loc,float *Latitude, float *Longitude)
  2388. /*  convert grid locator to latitude and longitude */
  2389. {
  2390.    char locator[6];
  2391.    float lat,lon;
  2392.    int i,nc;
  2393.  
  2394.    nc=strlen(loc);
  2395.  
  2396.    for ( i=0; i<=5; i++ ) {
  2397.       locator[i]=tolower(loc[i]);
  2398.    }
  2399.  
  2400.    for ( i=0; i<=1; i++ ) {
  2401.       if ((locator[i]<'a')||(locator[i]>'r')) return 1;
  2402.    }
  2403.  
  2404.    for ( i=2; i<=3; i++ ) {
  2405.       if ( (locator[i]<'0')||(locator[i]>'9') ) return 1;
  2406.    }
  2407.  
  2408.    if (nc > 4){
  2409.       for ( i=4; i<=5; i++ ) {
  2410.          if ( (locator[i]<'a')||(locator[i]>'x') ) return 1;
  2411.       }
  2412.    }
  2413.  
  2414.    lat = 10.0*(locator[1] - 'a') - 90.0;
  2415.    lat += locator[3] - '0';
  2416.  
  2417.    lon = 20.0*(locator[0] - 'a') - 180.0;
  2418.    lon += 2.0*(locator[2] - '0');
  2419.  
  2420.    if (nc > 4){
  2421.       lat += (locator[5] - 'a')/24.0 + 1.0/48.0;
  2422.       lon += (locator[4] - 'a')/12.0 + 1.0/24.0;
  2423.    }
  2424.  
  2425.    *Latitude  = lat;
  2426.    *Longitude = lon;
  2427.  
  2428.    return 0;
  2429. }
  2430. void Help(void)
  2431. {
  2432.    ClearTheScreen();
  2433.  
  2434.    printf (" You MUST specify an argument on command line for STP to do");
  2435.    printf (" anything\n");
  2436.    printf (" Enter STP and press ENTER to see list of commands available\n");
  2437.    printf ("\n");
  2438.    printf (" For DOS machines:\n");
  2439.    printf (" Make sure environment variable TZ is set in AUTOEXEC.BAT\n");
  2440.    printf (" Example for KD4QIO (Central Time Zone)\n");
  2441.    printf ("    SET TZ=CDT6CST\n");
  2442.    printf ("\n");
  2443.    printf (" On certain machines some of the calculations");
  2444.    printf (" requested take some time.\n");
  2445.    printf (" Please be patient....\n");
  2446.    printf ("\n");
  2447.    printf (" If all else fails try the following :\n");
  2448.    printf ("\n");
  2449.    printf (" Give me a call at (205) 837-5282 x1216 during the day\n");
  2450.    printf ("\n");
  2451.    printf (" mail:\n kd4qio@amsat.org\n harper@huntsville.sparta.com\n");
  2452.    printf ("\n or\n\n");
  2453.    printf (" C.Harper\n P.O. Box 18786\n Huntsville,Al 35804-8786\n");
  2454.    printf ("\n");
  2455. }
  2456. void HourMinuteSecond(float time,FILE *WhereToWrite,int DayOfTheWeek)
  2457. /*  convert time (in decimal hours) to day:hour:minute:second format */
  2458. {
  2459.    float Temptime;
  2460.  
  2461.    Temptime = time;
  2462.  
  2463.    TempDay  = ZERO;
  2464.    TempHour = AINT(Temptime);
  2465.  
  2466. /*  days */
  2467.  
  2468.    if (TempHour > 24.0){
  2469.       TempDay  = AINT(TempHour/24.0);
  2470.       TempHour = TempHour - 24.0*TempDay;
  2471.    }
  2472.  
  2473. /*  minutes */
  2474.    Temptime -= (24.0*TempDay+TempHour);
  2475.    Temptime *= 60.0;
  2476.  
  2477.    TempMin = AINT(Temptime);
  2478.   
  2479. /*  seconds */
  2480.    Temptime -= TempMin;
  2481.    Temptime *= 60.0;
  2482.  
  2483.    TempSec = AINT(Temptime);
  2484.  
  2485. /*  if DayOfTheWeek == -999 dont print anything */
  2486.  
  2487.    if (DayOfTheWeek == -999)return;
  2488.  
  2489. /* 
  2490.    set day of the week
  2491.  
  2492.    if DayOfTheWeek is < 0 then user does NOT want this info
  2493. */
  2494.    if (DayOfTheWeek >= 0){
  2495.       if (TempDay > ZERO){
  2496.          DayOfTheWeek += (int)TempDay;
  2497.          if (DayOfTheWeek > 6){
  2498.             DayOfTheWeek -= 7;
  2499.          }
  2500.       }
  2501.       fprintf (WhereToWrite,"%12s ",
  2502.       DayNames[DayOfTheWeek]);
  2503.    }else{
  2504.       fprintf (WhereToWrite,"%02.0f:",
  2505.       TempDay);
  2506.    }
  2507.  
  2508.    fprintf (WhereToWrite,"%02.0f:%02.0f:%02.0f",
  2509.       TempHour,
  2510.       TempMin,
  2511.       TempSec);
  2512. }
  2513. void InitConstants(void)
  2514. /* initialize constants */
  2515. {
  2516.  
  2517.    TRUE   = 1;
  2518.  
  2519.    ONE    = 1.0;
  2520.    ZERO   = 0.0;
  2521.    PI     = 3.141592654;
  2522.    TWOPI  = 2.0*PI;
  2523.    PIOVR2 = PI/2.0;
  2524.    D2R    = 180.0/PI;
  2525.    THIRD  = 1.0/3.0;
  2526.  
  2527.    RE     = 6378.140;       /* radius of the earth */
  2528.    GM     = 3.986E5;        /* earths gravitational constant km^3/s^2 */
  2529.    J2     = 1.08263E-3;     /* 2nd zonal coeff, earths gravity field */
  2530.  
  2531.    YM     = 365.25;         /* mean year, days */
  2532.    YT     = 365.2421938;    /* tropical year, days */
  2533.  
  2534.    YG     = 1992.0;     /*  year sun data from */
  2535.    G0     = 98.9205;
  2536.  
  2537.    MAS0 = 356.1230621;  /* mean anomaly for sun, deg */
  2538.    MASD = 0.98560027;   /* mean anomaly rate for sun, deg/day */
  2539.    INS  = 23.4403;      /* suns inclination, deg */
  2540.  
  2541.    EQC1 = 0.0334263;    /* suns equation of center terms */
  2542.    EQC2 = 0.00034915;   /* suns equation of center terms */
  2543.    EQC3 = 5.0e-6;       /* suns equation of center terms */
  2544.  
  2545.    AstroUnit = 1.495985e8; /*  one astronomical unit km */
  2546.  
  2547. /* 
  2548.    moon constants......
  2549.  
  2550.    Epoch.....
  2551. */
  2552.  
  2553.    Month=1.0;
  2554.    Day=0.0;
  2555.    Year=1993.0;
  2556.  
  2557.    MoonEpoch = DayNumber(Year,Month,Day);
  2558.  
  2559.    LO1 = 6.288750;
  2560.    LO2 =-1.274018;
  2561.    LO3 = 0.658309;
  2562.    LO4 = 0.213616;
  2563.    LO5 =-0.185596;
  2564.    LO6 =-0.114336;
  2565.    LO7 =-0.058793;
  2566.    LO8 =-0.057212;
  2567.    LO9 = 0.053320;
  2568.    LO10=-0.045874;
  2569.    LO11= 0.041024;
  2570.    LO12=-0.034718;
  2571.  
  2572.    LA1= 5.128189;
  2573.    LA2= 0.280606;
  2574.    LA3=-0.277693;
  2575.    LA4=-0.173238;
  2576.    LA5= 0.055413;
  2577.    LA6= 0.046272;
  2578.    LA7= 0.032573;
  2579.  
  2580.    HP0=0.950724;
  2581.    HP1=0.051818;
  2582.    HP2=0.009531;
  2583.    HP3=0.007843;
  2584.    HP4=0.002824;
  2585.    HP5=0.000857;
  2586.  
  2587.    LM0=359.6821;
  2588.    LMD=13.17640;
  2589.  
  2590.    MA0=201.2434;
  2591.    MAD=13.06499;
  2592.  
  2593.    F0 =99.2094;
  2594.    FD =13.22935;
  2595.    D0 =80.0093;
  2596.    DD =12.19075;
  2597.  
  2598.    RevolutionOld = 0.0;
  2599.    MonthOld      = 0;
  2600.    DayOld        = 0;
  2601.    HourOld       = 0.0;
  2602.    MinuteOld     = 0.0;
  2603.    SecondOld     = 0.0;
  2604.  
  2605. /*  length of line....how many dashes to print out in prediction mode */
  2606.  
  2607.    LengthOfLine = 80 - strlen(BLANKS) - 7;
  2608.  
  2609. /*  open file and read some control flags */
  2610.  
  2611.    if ((ControlFile = fopen (FlagFile,"r")) == 0){
  2612.       printf (" Cant open %s\n",FlagFile);
  2613.       exit(-1);
  2614.    }
  2615.  
  2616.    fgets (line1,80,ControlFile);
  2617.    sscanf (line1,"%f",&FindStep);
  2618.    FindStep /= 1440.0;
  2619.  
  2620.    fgets (line1,80,ControlFile);
  2621.    sscanf (line1,"%f",&FindAngleError);
  2622.    FindAngleError /= D2R;
  2623.  
  2624.    fgets (line1,80,ControlFile);
  2625.    sscanf (line1,"%f",&FindStepMinimum);
  2626.    FindStepMinimum /= 1440.0;
  2627.  
  2628.    fgets (line1,80,ControlFile);
  2629.    sscanf (line1,"%d",&PlotFlag);
  2630.  
  2631.    fgets (line1,80,ControlFile);
  2632.    sscanf (line1,"%d",&PlusMinus);
  2633.  
  2634.    fgets (line1,80,ControlFile);
  2635.    sscanf (line1,"%d",&GroundTrace);
  2636.  
  2637.    fgets (line1,80,ControlFile);
  2638.    sscanf (line1,"%f %f",&FrequencyModeB,&FrequencyModeS);
  2639.  
  2640.    fgets (line1,80,ControlFile);
  2641.    sscanf (line1,"%f %f",&BandWidthModeB,&BandWidthModeS);
  2642.  
  2643.    fgets (line1,80,ControlFile);
  2644.    sscanf (line1,"%f %f",&ReceiveGainModeB,&ReceiveGainModeS);
  2645.  
  2646.    fgets (line1,80,ControlFile);
  2647.    sscanf (line1,"%f %f",&ReceiveNoiseModeB,&ReceiveNoiseModeS);
  2648.  
  2649.    fgets (line1,80,ControlFile);
  2650.    sscanf (line1,"%f",&AverageNumberOfUsers);
  2651.  
  2652. /*  convert to dB */
  2653.  
  2654.    AverageNumberOfUsers = 10.0*log10(AverageNumberOfUsers);
  2655.  
  2656.    fclose (ControlFile);
  2657.  
  2658.    NumberOfPlots= - 1;
  2659.  
  2660. }
  2661. void InitAntenna(void)
  2662. /* initailize antenna unit vector in orbit plane coordinates */
  2663. {
  2664.    float COa,SOa,CLa,SLa;
  2665.  
  2666.    COa = cos(SATS[WhichSatellite].Alon);
  2667.    SOa = sin(SATS[WhichSatellite].Alon);
  2668.    CLa = cos(SATS[WhichSatellite].Alat);
  2669.    SLa = sin(SATS[WhichSatellite].Alat);
  2670.  
  2671.    SATS[WhichSatellite].ax = -CLa*COa;
  2672.    SATS[WhichSatellite].ay = -CLa*SOa;
  2673.    SATS[WhichSatellite].az = -SLa;
  2674. }
  2675. void InitQth(void)
  2676. /* process input Qth information */
  2677. {
  2678.  
  2679.    float D,SL,CL,Rx,Rz,SO,CO;
  2680.    int i;
  2681.  
  2682.    for (i=0;i<=NumberOfQth;i++){
  2683.       Qth[i].LA /=D2R;
  2684.       Qth[i].LO /=D2R;
  2685.       Qth[i].HT /=1000.;
  2686.  
  2687.       CL = cos(Qth[i].LA);
  2688.       SL = sin(Qth[i].LA);
  2689.       CO = cos(Qth[i].LO);
  2690.       SO = sin(Qth[i].LO);
  2691.  
  2692.       D = sqrt(RE2*CL*CL + ZZ*SL*SL);
  2693.       Rx = RE2/D + Qth[i].HT;
  2694.       Rz = ZZ/D + Qth[i].HT;
  2695.  
  2696. /* observers unit vectors UP EAST and NORTH in geocentric coords */
  2697.  
  2698.       Qth[i].Ux = CL*CO;
  2699.       Qth[i].Ex = -SO;
  2700.       Qth[i].Nx = -SL*CO;
  2701.       Qth[i].Uy = CL*SO;
  2702.       Qth[i].Ey = CO;
  2703.       Qth[i].Ny = -SL*SO;
  2704.       Qth[i].Uz = SL;
  2705.       Qth[i].Ez = ZERO;
  2706.       Qth[i].Nz = CL;
  2707.  
  2708. /* observers XYZ coords at earths surface */
  2709.  
  2710.       Qth[i].Ox = Rx*Qth[i].Ux;
  2711.       Qth[i].Oy = Rx*Qth[i].Uy;
  2712.       Qth[i].Oz = Rz*Qth[i].Uz;
  2713.  
  2714. /* observers velocity in geocentric coords (vz=0) */
  2715.  
  2716.       Qth[i].V0x = -Qth[i].Oy*W0;
  2717.       Qth[i].V0y =  Qth[i].Ox*W0;
  2718.       Qth[i].V0z = ZERO;
  2719.    }
  2720. }
  2721. void InitSatellite(void)
  2722. /* compute parameters associated with keplerian data */
  2723. {
  2724.    int i;
  2725.  
  2726.    for (i=0;i<=NumberOfSats;i++){
  2727.  
  2728.       RiseSet[i].RiseDay = ZERO;
  2729.       RiseSet[i].RiseHour = ZERO;
  2730.       RiseSet[i].RiseMinute = ZERO;
  2731.       RiseSet[i].RiseSecond = ZERO;
  2732.       RiseSet[i].SetDay = ZERO;
  2733.       RiseSet[i].SetHour = ZERO;
  2734.       RiseSet[i].SetMinute = ZERO;
  2735.       RiseSet[i].SetSecond = ZERO;
  2736.  
  2737.       SATS[i].RA /=D2R;
  2738.       SATS[i].IN /=D2R;
  2739.       SATS[i].WP /=D2R;
  2740.       SATS[i].MA /=D2R;
  2741.       SATS[i].MM *=TWOPI;
  2742.       SATS[i].M2 *=TWOPI;
  2743.  
  2744. /* convert satellite epoch to day number and fraction of a day */
  2745.  
  2746.       SATS[i].DE = DayNumber(SATS[i].YE,ONE,ZERO) + AINT(SATS[i].TE);
  2747.       SATS[i].TE = SATS[i].TE - AINT(SATS[i].TE);
  2748.  
  2749. /* mean motion, rad/s */
  2750.       SATS[i].N0 = SATS[i].MM/86400.0;            
  2751. /* semi major axis, km */
  2752.       SATS[i].A0 = (GM/SATS[i].N0/SATS[i].N0);
  2753. /*  
  2754.    if the sun has been input then force the semi-major axis to the 
  2755.    correct value.
  2756.    determine approximate period (minutes)
  2757. */
  2758.       if (strncmpi(SATS[i].SatsName,"sun",3) == 0){
  2759.          SATS[i].A0 = 149597870.0;
  2760.          SATS[i].Period = 365.0*24.0*60.0;
  2761.       }else{
  2762.          SATS[i].A0 = pow(SATS[i].A0,THIRD);         
  2763.          SATS[i].Period = 1.7374e-4*pow(SATS[i].A0,1.5);
  2764.       }
  2765.  
  2766. /*  escobal stuff 
  2767.       printf (" mean motion...file/1440  %15.8f",SATS[i].MM/1440.0);
  2768.       printf (" mean motion...escobal    %15.8f\n",
  2769.          0.07436574/pow(SATS[i].A0/RE,3.0/2.0));
  2770. */
  2771.  
  2772. /* semi minor axis, km */
  2773.       SATS[i].B0 = SATS[i].A0*sqrt(ONE-SATS[i].EC*SATS[i].EC);
  2774.  
  2775. /* sin & cos of inclination */
  2776.       SATS[i].SI = sin(SATS[i].IN);
  2777.       SATS[i].CI = cos(SATS[i].IN);
  2778.  
  2779. /* precession rate, rad/day */
  2780.       SATS[i].PC = RE*SATS[i].A0/(SATS[i].B0*SATS[i].B0);
  2781.       SATS[i].PC = 1.5*J2*SATS[i].PC*SATS[i].PC*SATS[i].MM;        
  2782.  
  2783. /* node precession rate, rad/day */
  2784.       SATS[i].QD = -SATS[i].PC*SATS[i].CI;                 
  2785.  
  2786. /* perigee precession rate, rad/day */
  2787.       SATS[i].WD = SATS[i].PC*(5.0*SATS[i].CI*SATS[i].CI-ONE)/2.0; 
  2788.  
  2789. /* drag coefficient */
  2790.       SATS[i].DC = -2.0*SATS[i].M2/SATS[i].MM/3.0;         
  2791.    }
  2792. }
  2793. void InitSun(void)
  2794. /*  bring sun data to current satellite epoch */
  2795. {
  2796.    float TEG;
  2797.  
  2798.    TEG = SATS[WhichSatellite].DE
  2799.         -DayNumber(YG,ONE,ZERO)
  2800.         +SATS[WhichSatellite].TE;
  2801.  
  2802.    GHAE = G0 + TEG*WE ;
  2803.    MRSE = G0 + TEG*WW + PI;
  2804.    MASE = (MAS0 + MASD*TEG)/D2R;
  2805. }
  2806. int InputCheck(int minans,int maxans)
  2807. /*  read input and check for answer between minans and maxans */
  2808. {
  2809.    int retvar;
  2810. loop:
  2811.    scanf("%d",&retvar);
  2812.    if (retvar < minans || retvar > maxans) {
  2813.       printf (" Improper Response!!!! %d\n",retvar);
  2814.       printf (" Should be between %d and %d\n",minans,maxans);
  2815.       printf (" Try Again : ");
  2816.       goto loop;
  2817.    }
  2818.    return retvar;
  2819. }
  2820. int KeyEvent(void)
  2821. {
  2822. #ifdef MSDOS
  2823.  
  2824.    int key = bioskey(1);
  2825.  
  2826.    if (key) key = bioskey(0);
  2827.  
  2828.    return key;
  2829.  
  2830. #else
  2831.  
  2832.    return bioskey();
  2833.  
  2834. #endif
  2835. }
  2836. float LinInt(int npt,float x[],float y[],float xin)
  2837. /*  linear interpolation routine */
  2838. {
  2839.    int i,ib;
  2840.    float scale,add;
  2841.  
  2842. /*  less than 1st value...return 1st */
  2843.    if (xin <= x[0]){
  2844.       return y[0];
  2845.    }
  2846.  
  2847. /*  greater than last value...return last */
  2848.    if (xin >= x[npt]){
  2849.       return y[npt];
  2850.    }
  2851.  
  2852. /*  interpolate.....linear */
  2853.  
  2854.    for (i=0;i<=npt;i++){
  2855.       if (xin <= x[i]){
  2856.          ib=i-1;
  2857.          scale=(xin-x[ib])/(x[i]-x[ib]);
  2858.          add=scale*(y[i]-y[ib]);
  2859.          return y[ib]+add;
  2860.       }
  2861.    }
  2862.  
  2863. /*  if failed for some reason....return 0.0 */
  2864.    return ZERO;
  2865. }
  2866. float LOSAngle(float AzimuthToCheck)
  2867. /*  determine minimum line of sight angle for home qth */
  2868. {
  2869.     float Angle;
  2870.    
  2871. /*  if only one angle then return that value */
  2872.  
  2873.    if (NumberOfObsAngles <= 0){
  2874.       return ObsAngle[0];
  2875.    }
  2876.  
  2877. /*  more than one value...interpolate */
  2878.    
  2879.    Angle=AzimuthToCheck;
  2880.    if (Angle < ZERO) Angle+=360.0;
  2881.  
  2882.    return LinInt(NumberOfObsAngles,ObsAzimuth,ObsAngle,Angle);
  2883. }
  2884. void MatchQth(char WhatToMatch[])
  2885. /*  read qth.dat and look for a match */
  2886. {
  2887.    FILE *OQthFile;
  2888.  
  2889.    char DummyLine[132],RestOfLine[80];
  2890.    char *p;
  2891.    float lat,lon,alt;
  2892.    int i,nc,Ndigit;
  2893.  
  2894.    if ((OQthFile = fopen(DxQthFile,"r")) == 0){
  2895.       printf (" could not open file %s\n",DxQthFile);
  2896.       exit(-1);
  2897.    } 
  2898.  
  2899.    for (i=0;i<strlen(WhatToMatch);i++){
  2900.       WhatToMatch[i]=tolower(WhatToMatch[i]);
  2901.    }
  2902.  
  2903.    while (fgets(DummyLine,132,OQthFile)){
  2904.  
  2905.       for (i=0;i<strlen(DummyLine);i++){
  2906.          DummyLine[i]=tolower(DummyLine[i]);
  2907.       }
  2908.  
  2909.       nc=strlen(DummyLine);
  2910.  
  2911. /*  handle cases where stophere and starthere are found */
  2912.  
  2913.       if(nc > 10){
  2914.  
  2915. /*  latitude */
  2916.          p=strtok(DummyLine," ");
  2917.          lat=atof(p);
  2918.  
  2919. /*  longitude */
  2920.          p=strtok(NULL," ");
  2921.          lon=atof(p);
  2922.  
  2923. /*  altitude */
  2924.          p=strtok(NULL," ");
  2925.          alt=atof(p);
  2926.  
  2927. /*  grid locator */
  2928. /*
  2929.          p=strtok(NULL," ");
  2930. */
  2931.  
  2932. /*  rest of string....*/
  2933.          p=strtok(NULL,"\n");
  2934.  
  2935.          if (strstr(p,WhatToMatch)){
  2936.  
  2937.             printf (" Matched %s ",p);
  2938.  
  2939.             if (strncmpi(GetStringFromUser(&Ndigit),"y",1) == 0){
  2940.                strncpy(line1,p,35);
  2941.                line1[35]='\0';
  2942.                AddQth(lat,lon,alt,line1);
  2943.             }
  2944.          }
  2945.       }
  2946.    } 
  2947.  
  2948.    fclose (OQthFile);
  2949. }
  2950. void MonthDay(float DN)
  2951. /*   
  2952.   return day of the week (monday,tuesday...) and
  2953.   month (january,febuary,...) given Day Number
  2954. */
  2955. {
  2956.    float D,Y,Dw,Mn;
  2957.  
  2958.    D   = DN;
  2959.    D  += 428.0;
  2960.    Dw  = fmod(D + 5.0, 7.0);
  2961.    Y   = AINT((D - 122.1)/YM);
  2962.    D  -= AINT(Y*YM);
  2963.    Mn  = AINT(D/30.61);
  2964.    D  -= AINT(Mn*30.61);
  2965.    Mn -= ONE;
  2966.  
  2967.    if (Mn > 12.0){
  2968.       Mn -= 12.0;
  2969.       Y  += ONE;
  2970.    }
  2971.  
  2972.    MonthPointer = (3*Mn-2)/3;
  2973.    DayPointer   = (int)Dw;
  2974.  
  2975.    LocalDay  = D;
  2976.    LocalYear = Y;
  2977. }
  2978. void MoonAzimuthElevation(int PrintOut)
  2979. /*  compute moon azimuth and elevation at current time */
  2980. {
  2981.    float T,LM,MA,MS,F,D,LS;
  2982.    float LOMoon,LAMoon,HPMoon;
  2983.    float XMoon,YMoon,ZMoon,Xm,Ym,Zm,DECMoon,RAMoon;
  2984.  
  2985.    float C,S;
  2986.    float xRx,xRy,xRz,xR,xU,xE,xN;
  2987.    float Hx,Hy,Hz;
  2988.  
  2989.    T  = DayNow - MoonEpoch + TimeNow;
  2990.  
  2991.    LM = (LM0 + LMD*T)/D2R;
  2992.    MA = (MA0 + MAD*T)/D2R;
  2993.    F  = (F0  + FD*T)/D2R;
  2994.    D  = (D0  + DD*T)/D2R;
  2995.  
  2996.    LOMoon  = LO1*sin(MA)+LO2*sin(MA-2.0*D)+LO3*sin(2.0*D)+LO4*sin(2.0*MA)
  2997.             +LO5*sin(MS)+LO6*sin(2.0*F)+LO7*sin(2.0*(MA-D))+LO8*sin(MA+MS-2.0*D)
  2998.             +LO9*sin(MA+2.0*D)+LO10*sin(MS-2.0*D)+LO11*sin(MA-MS)+LO12*sin(D);
  2999.  
  3000.    LOMoon = LM + LOMoon/D2R;
  3001.  
  3002.    LAMoon = LA1*sin(F)+LA2*sin(MA+F)+LA3*sin(F-MA)+LA4*sin(F-2.0*D)
  3003.            +LA5*sin(F+2.0*D-MA)+LA6*sin(2.0*D-MA-F)+LA7*sin(F+2.0*D); 
  3004.  
  3005.    LAMoon /= D2R;
  3006.  
  3007.    HPMoon = HP0+HP1*cos(MA)+HP2*cos(MA-2.0*D)+HP3*cos(2.0*D)
  3008.            +HP4*cos(2.0*MA)+HP5*cos(MA+2.0*D);
  3009.  
  3010.    HPMoon /= D2R;
  3011.  
  3012.    DistanceToMoon = (double)(RE/sin(HPMoon));
  3013.  
  3014.    XMoon = (float)DistanceToMoon*cos(LAMoon)*cos(LOMoon);
  3015.    YMoon = (float)DistanceToMoon*cos(LAMoon)*sin(LOMoon);
  3016.    ZMoon = (float)DistanceToMoon*sin(LAMoon);
  3017.  
  3018.    Xm = XMoon;
  3019.    Ym = YMoon*CNS - ZMoon*SNS;
  3020.    Zm = YMoon*SNS + ZMoon*CNS;
  3021.  
  3022.    C = cos(-GHAA);
  3023.    S = sin(-GHAA);
  3024.  
  3025. /*  moon vector in geocentric coordinates */
  3026.    Hx = Xm*C - Ym*S;
  3027.    Hy = Xm*S + Ym*C;
  3028.    Hz = Zm;
  3029.  
  3030. /*  azimuth/elevation */
  3031.  
  3032.    xRx = Hx-Qth[WhichObserver].Ox;
  3033.    xRy = Hy-Qth[WhichObserver].Oy;
  3034.    xRz = Hz-Qth[WhichObserver].Oz;
  3035.    xR = sqrt(xRx*xRx+xRy*xRy+xRz*xRz);
  3036.  
  3037.    xRx = xRx/xR;
  3038.    xRy = xRy/xR;
  3039.    xRz = xRz/xR;
  3040.  
  3041.    xU = xRx*Qth[WhichObserver].Ux
  3042.        +xRy*Qth[WhichObserver].Uy
  3043.        +xRz*Qth[WhichObserver].Uz;
  3044.  
  3045. /*  Qth[].Ez == ZERO */
  3046.  
  3047.    xE = xRx*Qth[WhichObserver].Ex
  3048.        +xRy*Qth[WhichObserver].Ey;
  3049.  
  3050.    xN = xRx*Qth[WhichObserver].Nx
  3051.        +xRy*Qth[WhichObserver].Ny
  3052.        +xRz*Qth[WhichObserver].Nz;
  3053.  
  3054.    CurrentMoonAz = ATAN2(xE,xN)*D2R;
  3055.    CurrentMoonEl = asin(xU)*D2R;
  3056.  
  3057.    if (PrintOut){
  3058.       printf (" Moon: %7.2f %7.2f %15.3f",
  3059.          CurrentMoonAz,
  3060.          CurrentMoonEl,
  3061.          DistanceToMoon);
  3062.    }
  3063. }
  3064. void MutualRealTime(int argv,char *argc[])
  3065. /*  determine mutual window real time */
  3066. {
  3067.    char UserInput[80];
  3068.    int i,nc;
  3069.    int KeyStroke;
  3070.    float CurrentDay,DayEnd,NumberOfDays,OutputRate;
  3071.    float lat,lon,alt;
  3072.  
  3073. /*  set outputfile to stdout for other routines */
  3074.  
  3075.    OutputFile = stdout;
  3076.  
  3077. /*  process qth parameters */
  3078.    InitQth();
  3079.  
  3080. /*  handle cases where user may not follow rules....
  3081.  
  3082.     2 arguments....need satellite name...force update rate to 1 second
  3083.     3 arguments....decide if name or update rate
  3084.     4 arguments....hope its all correct
  3085. */
  3086.  
  3087.    switch (argv){
  3088.       case 0:
  3089.       case 1:
  3090.       case 2:
  3091.          for (i=0;i<=NumberOfSats;i++){
  3092.             printf (" %s\n",SATS[i].SatsName);
  3093.          }
  3094.  
  3095.          printf (" Which Satellite from the list ?  ");
  3096.  
  3097.          gets(SatName);
  3098.  
  3099.          break;
  3100.  
  3101.       case 3:
  3102.          OutputRate=atof(argc[2]);
  3103.  
  3104. /*  
  3105.   if 3 arguments and numeric equivalent of argument is zero then
  3106.   assume it is a satellite name ... if not then assume user
  3107.   input update rate...get satellite name
  3108. */
  3109.          if (OutputRate == ZERO){
  3110.             strcpy(SatName,argc[2]);
  3111.             OutputRate = ONE;
  3112.          }else{
  3113.             for (i=0;i<=NumberOfSats;i++){
  3114.                printf (" %s\n",SATS[i].SatsName);
  3115.             }
  3116.  
  3117.             printf (" Which Satellite from the list ?  ");
  3118.  
  3119.             gets(SatName);
  3120.             nc=strlen(SatName);
  3121.  
  3122.          }
  3123.  
  3124.          break;
  3125.  
  3126.       case 4:
  3127.          strcpy(SatName,argc[2]);
  3128.          OutputRate = atof(argc[3]);
  3129.  
  3130.          break;
  3131.    } 
  3132.  
  3133.    nc=strlen(SatName);
  3134.  
  3135. /*  match satellite name with current data */
  3136.  
  3137.    WhichSatellite=-1;
  3138.    for (i=0;i<=NumberOfSats;i++){
  3139.       if (strnicmp(SatName,SATS[i].SatsName,nc) == 0){
  3140.          if (WhichSatellite == -1){
  3141.             WhichSatellite = i;
  3142.          }else{
  3143.             printf (" Multiple Names.  Please be more specific \n");
  3144.             exit(-1);
  3145.          }
  3146.       }
  3147.    }
  3148.  
  3149.    if (WhichSatellite < 0){
  3150.       printf (" Could Not Match Satellite %s in List\n",SatName);
  3151.       exit (-1);
  3152.    }else{
  3153.       printf (" Matched your input as %s\n",
  3154.          SATS[WhichSatellite].SatsName);
  3155.    }     
  3156.  
  3157. /*  initailize sun parameters */
  3158.  
  3159.       InitSun();
  3160.  
  3161. /*  initailize antenna parameters */
  3162.  
  3163.       InitAntenna();
  3164.  
  3165.    while (TRUE){
  3166. /*
  3167.   get local date for possible use as input
  3168. */
  3169.       UTCLOCTime();
  3170.  
  3171.       Month = (float)GMTmon + ONE;
  3172.       Day   = (float)GMTdaym;
  3173.       Year  = (float)GMTyear + BasicYear;
  3174.  
  3175.       Hour   = (float)GMThour;
  3176.       Minute = (float)GMTmin;
  3177.       Second = (float)GMTsec;
  3178.  
  3179. /*  look for keystroke */
  3180.       KeyStroke = KeyEvent();
  3181.  
  3182.       if (KeyStroke){
  3183. #ifdef MSDOS
  3184.          KeyStroke &= 0x00ff;
  3185. #endif
  3186. /*  user wants to enter an observer in the form of a grid locator */
  3187.  
  3188.          if ((char)KeyStroke == 'g'){
  3189.             ClearTheScreen();
  3190.  
  3191.             printf (" Grid Locator: ");
  3192.  
  3193.             strcpy(UserInput,GetStringFromUser(&nc));
  3194.  
  3195.             if (nc){
  3196.                if (GridLocator(UserInput,&lat,&lon) == 0){
  3197.                   AddQth(lat,lon,ZERO,UserInput);
  3198.                }else{
  3199.                   printf (" \a\aGrid Locator input error\n");
  3200.                }
  3201.             }
  3202.          }
  3203.  
  3204. /*  user wants to enter an observer in the form of lat,lon,alt */
  3205.          if ((char)KeyStroke == 'l'){
  3206.             ClearTheScreen();
  3207.  
  3208.             printf (" Latitude,Longitude,Altitude: ");
  3209.  
  3210.             strcpy(UserInput,GetStringFromUser(&nc));
  3211.             sscanf(UserInput,"%f %f %f",&lat,&lon,&alt);
  3212.             AddQth(lat,lon,alt,UserInput);
  3213.          }
  3214.  
  3215. /*  user wants to enter an observer name to match in qth.dat */
  3216.          if ((char)KeyStroke == 'o'){
  3217.             ClearTheScreen();
  3218.  
  3219.             printf (" City,Callsign or Grid: ");
  3220.  
  3221.             strcpy(UserInput,GetStringFromUser(&nc));
  3222.             MatchQth(UserInput);
  3223.          }
  3224.  
  3225. /* quit */
  3226.          if ((char)KeyStroke == 'q'){
  3227.             exit(0);
  3228.          }
  3229.       }
  3230.  
  3231. /*  get current day number */
  3232.       DayNow=DayNumber(Year,Month,Day);
  3233.  
  3234.       if (DeltTime(OutputRate)){
  3235.  
  3236.         ClearTheScreen();
  3237.  
  3238. /*  show user current local and GMT time */
  3239.  
  3240.          printf (" LOC %10s %10s %02d 19%02d %02d:%02d:%02d\n",
  3241.             DayNames[LOCdayw],
  3242.             MonthNames[LOCmon],
  3243.             LOCdaym,
  3244.             LOCyear,
  3245.             LOChour,
  3246.             LOCmin,
  3247.             LOCsec);
  3248.  
  3249.          printf (" GMT %10s %10s %02d 19%02d %02d:%02d:%02d\n",
  3250.             DayNames[GMTdayw],
  3251.             MonthNames[GMTmon],
  3252.             GMTdaym,
  3253.             GMTyear,
  3254.             GMThour,
  3255.             GMTmin,
  3256.             GMTsec);
  3257.  
  3258.          printf ("\n");
  3259.  
  3260. /*  current time in hours */
  3261.  
  3262.          TimeNow = (Hour+Minute/60.0+Second/3600.0)/24.0;
  3263.  
  3264. /*  get satellite location */
  3265.          SatellitePosition();
  3266.  
  3267. /*  get local information */
  3268.          WhichObserver = 0;
  3269.  
  3270. /*  determine azimuth and elevation from observers location */
  3271.          ObserverSatelliteGeometry();
  3272.  
  3273. /*  compute sun lit status */
  3274.  
  3275.          SunData();
  3276.  
  3277. /*  get mode info */
  3278.  
  3279.          WhatMode();
  3280.  
  3281.          printf (" Real Time Mutual for: %s",
  3282.             SATS[WhichSatellite].SatsName);
  3283.  
  3284.          printf (" MA: %5.1f Mode: %s",CurrentMAd,CurrentMode);
  3285.  
  3286. /*  sunlit status */
  3287.  
  3288.          printf (" Satel/SUN (%s)",ECL);
  3289.  
  3290.          printf ("\n");
  3291.          printf ("\n");
  3292.  
  3293.          printf ("  Az    El  Range Squnt  Doppler");
  3294.          printf ("  Status");
  3295.  
  3296.          printf ("\n");
  3297.  
  3298.          printf(" %4.0f %4.0f %6.0f",
  3299.             Azimuth*D2R,
  3300.             Elevation*D2R,
  3301.             RangeToSatellite);
  3302.  
  3303. /*  compute squint/rangerate */
  3304.  
  3305.          SquintRangeRate();
  3306.  
  3307. /*  if Omnidirectional antennas on ao-10/ao-13... then adjust squint angle */
  3308.  
  3309.          if (strncmpi(CurrentMode,"Bo",2) == 0) {
  3310.             SquintAngle = ABS(90.0 - SquintAngle);
  3311.          }
  3312.                    
  3313.          printf (" %4.0f",SquintAngle);
  3314.  
  3315. /*  compute doppler */
  3316.  
  3317.          Doppler(SATS[WhichSatellite].Beacon);
  3318.  
  3319.          printf ("   ");
  3320.  
  3321.          if (Elevation > ZERO){
  3322.             printf ("  Up   Home QTH");
  3323.          }else{
  3324.             printf (" Down  Home QTH");
  3325.          }
  3326.  
  3327.          printf ("\n");
  3328.  
  3329.          if (NumberOfQth > 0){
  3330.             for (i=1;i<=NumberOfQth;i++){
  3331.                WhichObserver = i;
  3332.  
  3333.                ObserverSatelliteGeometry();
  3334.  
  3335.                printf(" %4.0f %4.0f %6.0f",
  3336.                   Azimuth*D2R,
  3337.                   Elevation*D2R,
  3338.                   RangeToSatellite);
  3339.  
  3340. /*  compute squint/rangerate */
  3341.  
  3342.                SquintRangeRate();
  3343.  
  3344.                if (strncmpi(CurrentMode,"Bo",2) == 0) {
  3345.                   SquintAngle = ABS(90.0 - SquintAngle);
  3346.                }
  3347.  
  3348.                printf (" %4.0f",SquintAngle);
  3349.  
  3350. /*  compute doppler */
  3351.  
  3352.                Doppler(SATS[WhichSatellite].Beacon);
  3353.  
  3354.                printf ("   ");
  3355.  
  3356.                if (Elevation > ZERO){
  3357.                   printf ("  Up  ");
  3358.                }else{
  3359.                   printf (" Down ");
  3360.                }
  3361.  
  3362.                printf (" %s",
  3363.                   Qth[WhichObserver].QthLoc);
  3364.  
  3365.                printf ("\n");
  3366.             }
  3367.          }
  3368.  
  3369.          SecondOld = Second;
  3370.          MinuteOld = Minute;
  3371.          HourOld   = Hour;
  3372.       }
  3373.    }
  3374. }
  3375. void MutualViewing(int argv,char *argc[])
  3376. /*  determine mutual windows */
  3377. {
  3378.    int i,nc;
  3379.    float CurrentDay,DayEnd,NumberOfDays;
  3380.    float RangeLoc,AzimuthLoc,ElevationLoc,SquintLoc;
  3381.  
  3382. /*  read other qth data */
  3383.    ReadOtherQth();
  3384.  
  3385.    if (NumberOfQth == 0){
  3386.       printf (" NO Qth specified in qth.dat\n");
  3387.       exit(-1);
  3388.    }
  3389.  
  3390. /*  process qth parameters */
  3391.    InitQth();
  3392.  
  3393. /*  if user did not input satellite on command line then ask for input */
  3394.  
  3395.    if (argv <= 2){
  3396.       for (i=0;i<=NumberOfSats;i++){
  3397.          printf (" %s\n",SATS[i].SatsName);
  3398.       }
  3399.  
  3400.       printf (" Which Satellite from the list ?  ");
  3401.  
  3402.       gets(SatName);
  3403.    } else {
  3404.       strcpy(SatName,argc[2]);
  3405.    }
  3406.  
  3407.    nc=strlen(SatName);
  3408.  
  3409. /*  match satellite name with current data */
  3410.  
  3411.    WhichSatellite=-1;
  3412.    for (i=0;i<=NumberOfSats;i++){
  3413.       if (strnicmp(SatName,SATS[i].SatsName,nc) == 0){
  3414.          if (WhichSatellite == -1){
  3415.             WhichSatellite = i;
  3416.          }else{
  3417.             printf (" Multiple Names.  Please be more specific \n");
  3418.             exit(-1);
  3419.          }
  3420.       }
  3421.    }
  3422.  
  3423.    if (WhichSatellite < 0){
  3424.       printf (" Could Not Match Satellite %s in List\n",SatName);
  3425.       exit (-1);
  3426.    }else{
  3427.       printf (" Matched your input as %s\n",SATS[WhichSatellite].SatsName);
  3428.    }     
  3429.  
  3430. /*
  3431.   get local date for possible use as input
  3432. */
  3433.    UTCLOCTime();
  3434.  
  3435.    Month = (float)GMTmon + ONE;
  3436.    Day   = (float)GMTdaym;
  3437.    Year  = (float)GMTyear + BasicYear;
  3438.  
  3439. /*  get current day number */
  3440.    CurrentDay=DayNumber(Year,Month,Day);
  3441.  
  3442.    printf (" Enter Start Day (month day year) (<cr> for current date) : ");
  3443.  
  3444. /*  allow user option to hit return to get current date */
  3445.    strcpy(line1,GetStringFromUser(&nc));
  3446.  
  3447. /*  if nc > 0 then user input...otherwise use current information */
  3448.  
  3449.    if (nc > 0){
  3450.       sscanf(line1,"%f%f%f",&Month,&Day,&Year);
  3451.    }else{
  3452.       Month = (float)GMTmon + ONE;
  3453.       Day   = (float)GMTdaym;
  3454.       Year  = (float)GMTyear + BasicYear;
  3455.    }
  3456.  
  3457.    printf (" Number of Days for Prediction: ");
  3458.    NumberOfDays=GetFloatNumberFromUser(&nc);
  3459.  
  3460.    printf (" Step Size (minutes....input negative value for seconds): ");
  3461.    Step=GetFloatNumberFromUser(&nc);
  3462.  
  3463. /* 
  3464.    if Step is negative then user has input time step in seconds...
  3465.    convert to minutes 
  3466. */
  3467.  
  3468.    if (Step < ZERO){
  3469.       Step = ABS(Step)/60.0;
  3470.    }
  3471.  
  3472.    printf (" FileName to write output to (cr for screen): ");
  3473.  
  3474. /* read filename here */
  3475.  
  3476.    strcpy(FileName,GetStringFromUser(&nc));
  3477.  
  3478.    if (nc > 0){
  3479.       if ((OutputFile = fopen(FileName,"w")) == 0){
  3480.          printf ("Can't open file %s\n",FileName);
  3481.          exit(-1);
  3482.       }
  3483.    }else{
  3484.       OutputFile = stdout;
  3485.    }     
  3486.  
  3487.    DayStart = DayNumber(Year,Month,Day);
  3488.    DayEnd   = DayStart + NumberOfDays - ONE;
  3489.  
  3490. /*  initailize sun parameters */
  3491.  
  3492.    InitSun();
  3493.  
  3494. /*  initailize antenna parameters */
  3495.  
  3496.    InitAntenna();
  3497.  
  3498. /*  loop over time */
  3499.  
  3500.    RunningCounter    = 0;
  3501.    RunningCounterOld = 0;
  3502.  
  3503.    for (DayNow=DayStart;DayNow<=DayEnd;DayNow++){
  3504.       for (Hour=ZERO;Hour<=23.0;Hour++){
  3505.          for (Minute=ZERO;Minute<=60.0-Step;Minute=Minute+Step){
  3506.  
  3507.             RunningCounter++;
  3508.  
  3509.             TimeNow = (Hour+Minute/60.0)/24.0;
  3510.  
  3511. /*  get satellite location */
  3512.             SatellitePosition();
  3513.  
  3514. /*  get local information */
  3515.             WhichObserver = 0;
  3516.  
  3517. /*  determine azimuth and elevation from observers location */
  3518.             ObserverSatelliteGeometry();
  3519.  
  3520.             WhatMode();
  3521.  
  3522.             if (strncmpi(CurrentMode,"Bo",2) == 0) {
  3523.                SquintAngle = ABS(90.0 - SquintAngle);
  3524.             }
  3525.  
  3526. /* can we see it locally then check the other qth's */
  3527.  
  3528.             if (Elevation >= LOSAngle(Azimuth)){
  3529.  
  3530.                SunData();
  3531.                SquintRangeRate();
  3532.  
  3533.                RangeLoc     = RangeToSatellite;
  3534.                AzimuthLoc   = Azimuth;
  3535.                ElevationLoc = Elevation;
  3536.                SquintLoc    = SquintAngle;
  3537.  
  3538.                for (i=1;i<=NumberOfQth;i++){
  3539.  
  3540.                   WhichObserver = i;
  3541.                   ObserverSatelliteGeometry();
  3542.  
  3543.                   if (Elevation >= ZERO){
  3544.                      SunData();
  3545.                      SquintRangeRate();
  3546.  
  3547.                      PrintHeader(1,CurrentDay);
  3548.  
  3549.                      if (strncmpi(CurrentMode,"Bo",2) == 0) {
  3550.                         SquintAngle = ABS(90.0 - SquintAngle);
  3551.                         SquintLoc   = ABS(90.0 - SquintLoc);
  3552.                      }
  3553.  
  3554.                      fprintf (OutputFile,"%s",BLANKS);
  3555.  
  3556.                      if (Step < ONE){
  3557.                         Seconds = 60.0*(Minute-AINT(Minute));
  3558.                         Minutes = AINT(Minute);
  3559.                         fprintf (OutputFile,"%02.0f%02.0f%02.0f",
  3560.                            Hour,
  3561.                            Minutes,
  3562.                            Seconds);
  3563.                      }else{
  3564.                         fprintf (OutputFile,"%02.0f%02.0f",
  3565.                            Hour,
  3566.                            Minute);
  3567.                      }
  3568.  
  3569.                      fprintf (OutputFile," %4.0f %4.0f %7.0f %4.0f",
  3570.                         AzimuthLoc*D2R,
  3571.                         ElevationLoc*D2R,
  3572.                         RangeLoc,
  3573.                         SquintLoc);
  3574.  
  3575.                      fprintf (OutputFile,"  |%4.0f %4.0f %7.0f %4.0f",
  3576.                         Azimuth*D2R,
  3577.                         Elevation*D2R,
  3578.                         RangeToSatellite,
  3579.                         SquintAngle);
  3580.  
  3581.                      fprintf (OutputFile,"  %5.1f  %3s",
  3582.                         CurrentMAd,
  3583.                         CurrentMode);
  3584.  
  3585.                      fprintf (OutputFile," %s",
  3586.                         Qth[WhichObserver].QthLoc);
  3587.  
  3588.                      fprintf (OutputFile,"\n");
  3589.                   }
  3590.                }
  3591.             }
  3592.          }
  3593.  
  3594.          SecondOld = Second;
  3595.          MinuteOld = Minute;
  3596.          HourOld   = Hour;
  3597.       }
  3598.    }
  3599.  
  3600.    fclose (OutputFile);
  3601. }
  3602. void ObserverSatelliteGeometry(void)
  3603. /* determine azimuth and elevation from observer to satellite */
  3604. {
  3605.    Rx = Sx-Qth[WhichObserver].Ox;
  3606.    Ry = Sy-Qth[WhichObserver].Oy;
  3607.    Rz = Sz-Qth[WhichObserver].Oz;
  3608.  
  3609.    RangeToSatellite = sqrt(Rx*Rx+Ry*Ry+Rz*Rz);
  3610.  
  3611.    Rx = Rx/RangeToSatellite;
  3612.    Ry = Ry/RangeToSatellite;
  3613.    Rz = Rz/RangeToSatellite;
  3614.  
  3615.    U = Rx*Qth[WhichObserver].Ux
  3616.       +Ry*Qth[WhichObserver].Uy
  3617.       +Rz*Qth[WhichObserver].Uz;
  3618.    E = Rx*Qth[WhichObserver].Ex
  3619.       +Ry*Qth[WhichObserver].Ey;
  3620.    N = Rx*Qth[WhichObserver].Nx
  3621.       +Ry*Qth[WhichObserver].Ny
  3622.       +Rz*Qth[WhichObserver].Nz;
  3623.  
  3624.    Azimuth   = ATAN2(E,N);
  3625.    Elevation = asin(U);
  3626. }
  3627. float P0to360(float ans)
  3628. /*  Insure ans between 0 and 360 */
  3629. {
  3630.    float n360,m360;
  3631.  
  3632.    if (ans < 0.0){
  3633.       n360=ABS(ans)/360.0;
  3634.       m360=AINT(n360);
  3635.       return ans+m360*360.0;
  3636.    }   
  3637.    if (ans > 360.0){
  3638.       n360=ans/360.0;
  3639.       m360=AINT(n360);
  3640.       return ans-m360*360.0;
  3641.    }
  3642.    return ans;
  3643. }
  3644. void PercivedSignalStrength(FILE *WhereToWrite)
  3645. /*  estimate signal strength */
  3646. {
  3647.    float Gain,SofT,Freq,PathLoss,Ps,Pn,Tsky,SquintTerm;
  3648.  
  3649. /*  only perform this calculation for ao-13 and ao-10 */
  3650.  
  3651.    if (strncmpi(SATS[WhichSatellite].SatsName,"ao-13",5) != 0 && 
  3652.        strncmpi(SATS[WhichSatellite].SatsName,"oscar10",7) != 0 &&
  3653.        strncmpi(SATS[WhichSatellite].SatsName,"oscar 10",8) != 0 &&
  3654.        strncmpi(SATS[WhichSatellite].SatsName,"ao13sm",6) != 0){
  3655.  
  3656.       return;
  3657.    }
  3658.  
  3659. /*
  3660.    Methodology taken from Satellite Experimenters Handbook
  3661.  
  3662.    Predicting signal levels pg 13-12
  3663.  
  3664.    assumptions:
  3665.  
  3666.       antenna pattern modeled using 2(n+1)cos(theta)^n
  3667.  
  3668.        n   Gain    Half-Power
  3669.            (dBi)   BeamWidth
  3670.       -------------------------
  3671.        0   3.0       180
  3672.       0.5  4.8       151
  3673.        1   6.0       120
  3674.        2   7.8        90
  3675.        3   9.0        74.9
  3676.        4  10.0        65.5
  3677.        6  11.5        54.0
  3678.        8  12.6        47.0
  3679.       12  14.1        38.6
  3680.       16  15.3        33.5
  3681.       20  16.2        30.0
  3682.       24  17.0        27.4
  3683.  
  3684.  
  3685.    ao-13 parameters
  3686.  
  3687.    Mode B
  3688.      Gain is 6.0 dbi
  3689.      Half-Power Beamwidth 100 deg
  3690.      Power is 12.5 W average
  3691.      n = 1.5 (match beamwidth)
  3692.   
  3693.    Mode S
  3694.      Gain is 14.2 dbi
  3695.      Half-Power Beamwidth 39 deg
  3696.      Power is 1.25 W continous
  3697.      n = 12 (match beamwidth)
  3698. */
  3699.  
  3700. /*  set gain factor for mode B antennas */
  3701.  
  3702.    if (strncmpi(CurrentMode,"Bo",2) == 0) {
  3703.       Gain=0.1;
  3704.    } else {
  3705.       Gain=1.5;
  3706.    }
  3707.  
  3708. /*  limit squint angles to 90 degrees or less */
  3709.  
  3710.    if (SquintAngle < 89.99){
  3711.       SquintTerm = log10(cos(SquintAngle/D2R));
  3712.    }else{
  3713.       SquintTerm = -4.7;
  3714.    }
  3715.  
  3716.     
  3717. /*  if Mode B, Mode Bo or Mode BS then compute gain and write out */
  3718.  
  3719.    SofT = ZERO;
  3720.  
  3721.    if (strncmpi(CurrentMode,"B",1) == 0){
  3722.  
  3723. /*  path loss at mode b frequency */
  3724.  
  3725.       PathLoss = 32.4 
  3726.                 +20.0*log10(FrequencyModeB) 
  3727.                 +20.0*log10(RangeToSatellite);
  3728.  
  3729. /* 
  3730.    40.969  is 10.0*log10(12.5*1000.0) 
  3731.  
  3732.    33.010  is 10.0*log10(2000)  follow same logic as G3RUH for Mode S
  3733.  
  3734.    This assumes you are only person on satellite.....
  3735.  
  3736.    Correct for number of users by subtracting off...
  3737.       10.0*log10(Number Of Users)
  3738. */
  3739.  
  3740.       Ps = 33.010 - AverageNumberOfUsers
  3741.           +10*log(2.0*Gain+ONE)
  3742.           +10.0*Gain*SquintTerm
  3743.           +ReceiveGainModeB
  3744.           -PathLoss;
  3745.  
  3746. /*  -198.79  is 10.0*log10(boltzmans constant in mw) */
  3747.  
  3748.       Pn = -198.79
  3749.            +10.0*log10(ReceiveNoiseModeB + TemperatureOfSky(FrequencyModeB))
  3750.            +10.0*log10(BandWidthModeB);
  3751.  
  3752.       SofT = Ps - Pn;
  3753.       
  3754.       fprintf (WhereToWrite," %8.2f ",SofT);
  3755.    }
  3756.  
  3757.    if (PlotFlag)fprintf (PlotFile1,"\t%8.2f",SofT);
  3758.  
  3759. /*  perform Mode S calculation in case BS or S only  */
  3760.  
  3761.    PathLoss = 32.4 
  3762.              +20.0*log10(FrequencyModeS) 
  3763.              +20.0*log10(RangeToSatellite);
  3764.  
  3765. /* 
  3766.    30.969  is 10.0*log10(1.25*1000.0)
  3767.  
  3768.    23.010  is 10.0*log10(200)   G3RUH (10/1/94)... avg power
  3769.  
  3770.    This assumes you are only person on satellite.....
  3771.  
  3772.    Correct for number of users by subtracting off...
  3773.       10.0*log10(Number Of Users)
  3774.  
  3775.    13.9794 is 10.0*log (2*12+1)
  3776. */
  3777.  
  3778.    Ps = 23.010 - AverageNumberOfUsers
  3779.        +13.9794
  3780.        +120.0*SquintTerm
  3781.        +ReceiveGainModeS
  3782.        -PathLoss;
  3783.  
  3784. /*  -198.79  is 10.0*log10(boltzmans constant in mw) */
  3785.  
  3786.    Pn = -198.79
  3787.         +10.0*log10(ReceiveNoiseModeS + TemperatureOfSky(FrequencyModeS));
  3788.         +10.0*log10(BandWidthModeS);
  3789.  
  3790.    SofT = Ps - Pn;
  3791.  
  3792. /*  if Mode S only then format slightly different from Mode BS */
  3793.  
  3794.    if (strncmpi(CurrentMode,"S..",3) == 0){
  3795.       fprintf (WhereToWrite,"           %8.2f ",SofT);
  3796.    }
  3797.    if (strncmpi(CurrentMode,"BS.",3) == 0){
  3798.       fprintf (WhereToWrite," %8.2f ",SofT);
  3799.    }
  3800.  
  3801.    if (PlotFlag)fprintf (PlotFile1,"\t%8.2f\n",SofT);
  3802.  
  3803. }
  3804. void PrintHeader(int OptionFlag,float CurrentDay)
  3805. /*  print header */
  3806. {
  3807.    int TimeLag;
  3808.  
  3809. /*  
  3810.   convert daynumber to a month (January,February,...) and 
  3811.   day (Monday,Tuesday,...)
  3812. */
  3813.    MonthDay(DayNow);
  3814.  
  3815. /*  
  3816.   look for cases where we have had a gap in prediction
  3817.   so we can write a new header block
  3818. */
  3819.  
  3820.    TimeLag = 0;
  3821.  
  3822.    if (RunningCounter - RunningCounterOld > 1)TimeLag = 1;
  3823.  
  3824.    RunningCounterOld = RunningCounter;
  3825.  
  3826.    if (MonthPointer != MonthOld || DayPointer != DayOld || TimeLag){
  3827.  
  3828.       if (PlotFlag){
  3829.  
  3830.          NumberOfPlots++;
  3831.  
  3832. /*  close up old files */
  3833.  
  3834.          if (NumberOfPlots > 0){
  3835.             fclose (PlotFile1);
  3836.             fclose (PlotFile2);
  3837.          }
  3838.  
  3839. /*  azimuth,elevation.... data */
  3840.  
  3841.          sprintf(PlotTemp,"azel%d.dat",NumberOfPlots);
  3842.          PlotFile1 = fopen (PlotTemp,"w");
  3843.  
  3844.          fprintf (PlotFile1,"9\nTime\nHour/Minute\nAz\nEl\nLat\nLon\nAlt\n");
  3845.          fprintf (PlotFile1,"S/N B\nS/N S\n");
  3846.  
  3847. /*  sub satellite footprint */
  3848.  
  3849.          sprintf(PlotTemp,"sbst%d.dat",NumberOfPlots);
  3850.          PlotFile2 = fopen (PlotTemp,"w");
  3851.  
  3852.          fprintf (PlotFile2,"4\nLat\nLon\nAlt\nTime\n");
  3853.  
  3854. /*  rise/set data */
  3855.  
  3856.          sprintf(PlotTemp,"rsst%d.dat",NumberOfPlots);
  3857.          PlotFile3 = fopen (PlotTemp,"w");
  3858.  
  3859.          fprintf (PlotFile3,"6\nTime\nAz\nEl\nLat\nLon\nAlt\n");
  3860.  
  3861.          fprintf (PlotFile3,"%02.0f%02.0f %7.2f %7.2f %10.3f %10.3f",
  3862.             Hour,Minute,
  3863.             Azimuth*D2R,
  3864.             Elevation*D2R,
  3865.             ATAN2(Sy,Sx)*D2R,
  3866.             asin(Sz/AltitudeOfSatellite)*D2R);
  3867.  
  3868.          if (GroundTrace){
  3869.             fprintf (PlotFile3,"%10.3f\n",ZERO);
  3870.          }else{
  3871.             fprintf (PlotFile3,"%10.3f\n",AltitudeOfSatellite-RE);
  3872.          }
  3873.  
  3874.          fclose (PlotFile3);
  3875.       }
  3876.  
  3877.       RevolutionOld = RevolutionNumber;
  3878.       MonthOld      = MonthPointer;
  3879.       DayOld        = DayPointer;
  3880.  
  3881.       fprintf (OutputFile,"%s",BLANKS);
  3882.       DashLine(OutputFile);
  3883.  
  3884.       fprintf (OutputFile,"%s",BLANKS);
  3885.       fprintf (OutputFile,"Prediction Date : %s %s %02.0f %4.0f",
  3886.          DayNames[DayPointer],
  3887.          MonthNames[MonthPointer],
  3888.          LocalDay,
  3889.          LocalYear);
  3890.       fprintf (OutputFile,"    Amsat Day: %6.0f\n",
  3891.               DayNow - 722100.0);
  3892.  
  3893.       fprintf (OutputFile,"%s",BLANKS);
  3894.  
  3895. /*  Qth[].QthLoc from qth.dat file will not have a carriage return at end */
  3896.  
  3897.       switch (OptionFlag){
  3898.          case 0:
  3899.             fprintf (OutputFile,"Prediction For  : %s",
  3900.                Qth[WhichObserver].QthLoc);
  3901.             break;
  3902.  
  3903.          case 1:
  3904.             if (NumberOfQth == 1){
  3905.                fprintf (OutputFile,"Prediction For  : %s\n",
  3906.                   Qth[WhichObserver].QthLoc);
  3907.             }else{
  3908.                fprintf (OutputFile,"Prediction For  : %2d Sites\n",
  3909.                   NumberOfQth);
  3910.             }
  3911.             break;
  3912.       }
  3913.  
  3914.       fprintf (OutputFile,"%s",BLANKS);
  3915.       fprintf (OutputFile,"Performed on    : %s %02d 19%02d %02d:%02d:%02d\n",
  3916.          MonthNames[LOCmon],
  3917.          LOCdaym,
  3918.          LOCyear,
  3919.          LOChour,
  3920.          LOCmin,
  3921.          LOCsec);
  3922.  
  3923.       fprintf (OutputFile,"%s",BLANKS);
  3924.       fprintf (OutputFile,"Satellite       : %s",
  3925.          SATS[WhichSatellite].SatsName);
  3926.  
  3927.       fprintf (OutputFile," (Set:%6d, Rev #:%6.0f, Age:%6.1f Days)",
  3928.          SATS[WhichSatellite].EleSet,
  3929.          RevolutionNumber,
  3930.          CurrentDay-(SATS[WhichSatellite].DE+SATS[WhichSatellite].TE));
  3931.       fprintf (OutputFile,"\n");
  3932.  
  3933.       fprintf (OutputFile,"%s",BLANKS);
  3934.       fprintf (OutputFile,"Illumination    : %3.1f%%  SAZ/SEL %4.1f/%4.1f ",
  3935.          Illumination,
  3936.          SunAzimuth*D2R,
  3937.          SunElevation*D2R);
  3938.       fprintf (OutputFile,"AP/RAAN %4.1f/%4.1f \n",
  3939.          AP*D2R,
  3940.          RAAN*D2R);
  3941.  
  3942.       fprintf (OutputFile,"%s",BLANKS);
  3943.       fprintf (OutputFile,"Alon/Alat       : %4.1f/%4.1f  ",
  3944.          SATS[WhichSatellite].Alon*D2R,SATS[WhichSatellite].Alat*D2R);
  3945.       fprintf (OutputFile,"Frequency: %6.4f MHz \n",
  3946.          SATS[WhichSatellite].Beacon*1.0e-6);
  3947.  
  3948.       fprintf (OutputFile,"%s",BLANKS);
  3949.       DashLine(OutputFile);
  3950.  
  3951.       switch (OptionFlag){
  3952.          case 0:
  3953.             fprintf (OutputFile,"%s",BLANKS);
  3954.             fprintf (OutputFile," UTC");
  3955.             if (Step < ONE)fprintf (OutputFile,"  ");
  3956.             fprintf (OutputFile,"   Az   El   Range  Sqnt    ");
  3957.             fprintf (OutputFile,"Dopplr    MA  Mode");
  3958.             fprintf (OutputFile,"  Sun    PRSL");
  3959.             fprintf (OutputFile,"\n");
  3960.  
  3961.             fprintf (OutputFile,"%s",BLANKS);
  3962.             DashLine(OutputFile);
  3963.  
  3964.             break;
  3965.  
  3966.          case 1:
  3967.             fprintf (OutputFile,"%s",BLANKS);
  3968.             fprintf (OutputFile,"            Local");
  3969.             if (Step < ONE)fprintf (OutputFile,"  ");
  3970.             fprintf (OutputFile,"            |");
  3971.             fprintf (OutputFile,"             Dx");
  3972.             fprintf (OutputFile,"\n");
  3973.  
  3974.             fprintf (OutputFile,"%s",BLANKS);
  3975.             fprintf (OutputFile," UTC");
  3976.             if (Step < ONE)fprintf (OutputFile,"  ");
  3977.             fprintf (OutputFile,"   Az   El   Range  Sqnt");
  3978.             fprintf (OutputFile," |");
  3979.             fprintf (OutputFile,"  Az   El   Range  Sqnt");
  3980.             fprintf (OutputFile,"   MA   Mde");
  3981.             fprintf (OutputFile,"\n");
  3982.  
  3983.             fprintf (OutputFile,"%s",BLANKS);
  3984.             DashLine(OutputFile);
  3985.             break;
  3986.  
  3987.          default:
  3988.             printf (" Cant match option %d in PrintHeader\n",
  3989.                OptionFlag);
  3990.             exit(-1);
  3991.       }
  3992.    }
  3993. }
  3994. void ReadKepler(void)
  3995. /*  read kepler data in kepler text format */
  3996. {
  3997.    char line[80],string1[80],string2[80],string3[80];
  3998.  
  3999.    int i,imat,nc,EpochYear;
  4000.    long int SatNum,EleSet,EpochRev;
  4001.  
  4002.    float EpochDay,Inclination,RAofAN,Eccentricity,ArgPerigee;
  4003.    float MeanMotion,DecayRate;
  4004.  
  4005. /*  get lengths of key words */
  4006.  
  4007.    for (i=0;i<=MaxKeyWords;i++){
  4008.       LenKeyWords[i]=strlen(KeyWords[i]);
  4009.    }
  4010.  
  4011.    NumberOfSats=-1;
  4012.  
  4013. /*  
  4014.    read through file.
  4015.    try to match to one of keywords.
  4016. */
  4017.  
  4018.    while (fgets(line,80,KeplerDataFile)){
  4019.       imat=0;
  4020.  
  4021.       for (i=0;i<=MaxKeyWords;i++){
  4022.          if (strncmpi(line,KeyWords[i],LenKeyWords[i]) == 0){
  4023.             imat=1;
  4024.             switch (i){
  4025.                case 0:
  4026.                   NumberOfSats++;
  4027.  
  4028.                   if (NumberOfSats > MaxSat) {
  4029.                      NumberOfSats--;
  4030.                      printf (" Read in Maximum Number of Satellites %d\n",
  4031.                                     MaxSat);
  4032.                      printf (" Last Satellite Read in is %s\n",
  4033.                                     SATS[NumberOfSats].SatsName);
  4034.                      return;
  4035.                   }
  4036.  
  4037.                   sscanf(line,"%s %s",&string1,&string2);
  4038.  
  4039.                   nc=RealEndOfLine(string2);
  4040.                   strncpy(SATS[NumberOfSats].SatsName,string2,nc);
  4041. /*
  4042.                   SATS[NumberOfSats].SatsName = strdup(string2);
  4043. */
  4044.  
  4045. /*  set default values for mode information */
  4046.  
  4047.                   SATS[NumberOfSats].SquintOpt = 0;
  4048.                   SATS[NumberOfSats].NumModes = 0;
  4049.  
  4050.                   Modes[NumberOfSats][0].MinPhase = 0;
  4051.                   Modes[NumberOfSats][0].MaxPhase = 256;
  4052.  
  4053.                   strcpy(Modes[NumberOfSats][0].ModeStr,"---");
  4054.  
  4055. /*  set default values for beacon frequency, alat and alon */
  4056.  
  4057.                   SATS[NumberOfSats].Beacon = 145.00 * 1.0e6;
  4058.                   SATS[NumberOfSats].Alat   = 0.0;
  4059.                   SATS[NumberOfSats].Alon   = 0.0;
  4060.  
  4061.                   break;
  4062.      
  4063.                case 1:
  4064.                   sscanf(line,"%s %s %ld",&string1,&string2,&SatNum);
  4065.                   SATS[NumberOfSats].SatNum=SatNum;
  4066.                   break;
  4067.  
  4068.                case 2:
  4069.                   sscanf(line,"%s %s %2d %f",
  4070.                      &string1,&string2,&EpochYear,&EpochDay);
  4071.                   SATS[NumberOfSats].YE = BasicYear + EpochYear;
  4072.                   SATS[NumberOfSats].TE = EpochDay;
  4073.                   break;
  4074.  
  4075.                case 3:
  4076.                   sscanf(line,"%s %s %ld",&string1,&string2,&EleSet);
  4077.                   SATS[NumberOfSats].EleSet=EleSet;
  4078.                   break;
  4079.  
  4080.                case 4:
  4081.                   sscanf(line,"%s %f",&string1,&Inclination);
  4082.                   SATS[NumberOfSats].IN = Inclination;
  4083.                   break;
  4084.  
  4085.                case 5:
  4086.                   sscanf(line,"%s %s %s %f",
  4087.                      &string1,&string2,&string3,&RAofAN);
  4088.                   SATS[NumberOfSats].RA = RAofAN;
  4089.                   break;
  4090.  
  4091.                case 6:
  4092.                   sscanf(line,"%s %f",&string1,&Eccentricity);
  4093.                   SATS[NumberOfSats].EC = Eccentricity;
  4094.                   break;
  4095.  
  4096.                case 7:
  4097.                   sscanf(line,"%s %s %s %f",
  4098.                      &string1,&string2,&string3,&ArgPerigee);
  4099.                   SATS[NumberOfSats].WP = ArgPerigee;
  4100.                   break;
  4101.  
  4102.                case 8:
  4103.                   sscanf(line,"%s %s %f",&string1,&string2,&MeanAnomaly);
  4104.                   SATS[NumberOfSats].MA = MeanAnomaly;
  4105.                   break;
  4106.  
  4107.                case 9:
  4108.                   sscanf(line,"%s %s %f",&string1,&string2,&MeanMotion);
  4109.                   SATS[NumberOfSats].MM = MeanMotion;
  4110.                   break;
  4111.  
  4112.                case 10:
  4113.                   sscanf(line,"%s %s %f",&string1,&string2,&DecayRate);
  4114.                   SATS[NumberOfSats].M2 = DecayRate;
  4115.                   break;
  4116.  
  4117.                case 11:
  4118.                   sscanf(line,"%s %s %ld",&string1,&string2,&EpochRev);
  4119.                   SATS[NumberOfSats].RV = EpochRev;
  4120.                   break;
  4121.  
  4122.                case 12:
  4123.                   break;
  4124.             }
  4125.          }
  4126.       }
  4127.       if (imat == 0 && strlen(line) > 1){
  4128.          printf (" Could Not Match KeyWord %s",line);
  4129.          exit(-1);
  4130.       }
  4131.    }
  4132. }
  4133. void ReadKeps(void)
  4134. /*  
  4135.   control routine to read in keplerian data
  4136. */
  4137. {
  4138. /*  open data file
  4139.     nasa.dat is searched for 1st
  4140.     kepler.dat is next
  4141. */
  4142.  
  4143.    if ((NasaDataFile = fopen(NasaFile,"r")) == 0){
  4144.       if ((KeplerDataFile = fopen(KeplerFile,"r")) == 0){
  4145.          printf (" Cant open %s or %s !!!\n",
  4146.             NasaFile,KeplerFile);
  4147.          exit(-1);
  4148.       }else{
  4149.          ReadKepler();
  4150.       }
  4151.    }else{
  4152.       ReadNasa();
  4153.    }
  4154.  
  4155. /*  correct revolution counter */
  4156.  
  4157.    RevolutionCorrection();
  4158.  
  4159. /*  initialize satellite values */
  4160.  
  4161.    InitSatellite();
  4162. }
  4163. void ReadLocalQth(void)
  4164. {
  4165. /*  read in Qth information */
  4166.  
  4167.    FILE *QthFile;
  4168.  
  4169.    char line[80];
  4170.  
  4171.    int i,n;
  4172.    float adum,bdum;
  4173.  
  4174.    NumberOfQth = 0;
  4175.  
  4176. /*  if file does not exists then get info and create it */
  4177.  
  4178.    if ((QthFile = fopen(MyQthFile,"r")) == 0){
  4179.       printf (" STP initialization....Qth information\n\n");
  4180.  
  4181.       printf (" Your Latitude (Degrees): ");
  4182.       Qth[0].LA=GetFloatNumberFromUser(&n);
  4183.  
  4184.       printf (" Your Longitude (Degrees): ");
  4185.       Qth[0].LO=GetFloatNumberFromUser(&n);
  4186.  
  4187.       printf (" Your Altitude (Meters): ");
  4188.       Qth[0].HT=GetFloatNumberFromUser(&n);
  4189.  
  4190.       printf (" Identifier (Callsign, etc): ");
  4191.       strcpy(Qth[0].QthLoc,GetStringFromUser(&n));
  4192.  
  4193.       printf (" # of entries in minimum look angle table (at least 1): ");
  4194.  
  4195.       NumberOfObsAngles=InputCheck(1,MaxObsAng);
  4196.  
  4197.       printf (" input values seperated by a space \n");
  4198.  
  4199.       for (i=0;i<NumberOfObsAngles;i++){
  4200.          printf (" azimuth # %3d , minimum angle ",i+1);
  4201.          scanf("%f%f",&ObsAzimuth[i],&ObsAngle[i]);
  4202.          gets(line);
  4203.       }
  4204.  
  4205.       QthFile = fopen(MyQthFile,"a");
  4206.  
  4207.       fprintf (QthFile,"%s\n",Qth[0].QthLoc);
  4208.       fprintf (QthFile,"%12.6f\n",Qth[0].LA);
  4209.       fprintf (QthFile,"%12.6f\n",Qth[0].LO);
  4210.       fprintf (QthFile,"%12.6f\n",Qth[0].HT);
  4211.  
  4212.       for (i=0;i<NumberOfObsAngles;i++){
  4213.          fprintf (QthFile,"%10.2f %10.2f\n",ObsAzimuth[i],ObsAngle[i]);
  4214.       }
  4215.  
  4216.       fclose (QthFile);
  4217.  
  4218.       return;
  4219.    }
  4220.  
  4221. /*  file exists....read in data */
  4222.  
  4223.    fgets (Qth[0].QthLoc,40,QthFile);
  4224.  
  4225.    fgets (line,80,QthFile);
  4226.    sscanf(line,"%f",&Qth[0].LA);
  4227.  
  4228.    fgets (line,80,QthFile);
  4229.    sscanf(line,"%f",&Qth[0].LO);
  4230.  
  4231.    fgets (line,80,QthFile);
  4232.    sscanf(line,"%f",&Qth[0].HT);
  4233.  
  4234. /*  read in minimum line of sight angle */
  4235.  
  4236.    NumberOfObsAngles=-1;
  4237.  
  4238.    while(fgets(line,80,QthFile)){
  4239.       sscanf(line,"%f %f",&adum,&bdum);
  4240.  
  4241.       if(adum < ZERO)adum=360.0+adum;
  4242.  
  4243.       NumberOfObsAngles++;
  4244.  
  4245. /*  convert to radians */
  4246.  
  4247.       ObsAzimuth[NumberOfObsAngles] = adum/D2R;
  4248.       ObsAngle[NumberOfObsAngles]   = bdum/D2R;
  4249.  
  4250.       if (bdum < ZERO)return;
  4251.  
  4252.       if (NumberOfObsAngles > MaxObsAng){
  4253.          printf (" Too many angles in Qth file !!!!\n");
  4254.          exit (-1);
  4255.       }
  4256.    }
  4257. }
  4258. void ReadMode(void)
  4259. /*  read satellite mode file */
  4260. {
  4261.    char dum[11],mode[3];
  4262.    int i,min,max,imatch,whichsat,nc;
  4263.    float frqbec;
  4264.  
  4265.    if ((ModeDataFile = fopen(ModeFile,"r")) == 0){
  4266.       printf(" %s not found......use default data\n",ModeFile);
  4267.       return;
  4268.    }
  4269.  
  4270. /*
  4271.     read file and look for record... satellite:
  4272.     when found check for match to current sat
  4273. */
  4274.  
  4275.    imatch = 0;
  4276.    whichsat = -1;
  4277.  
  4278.    while (fgets(line1,80,ModeDataFile)){
  4279.  
  4280.       if(strnicmp(line1,"satellite:",10) == 0){
  4281.          imatch = 0;
  4282.      sscanf(line1,"%s%s",dum,SatName);
  4283.  
  4284. /*  try to match satellite name */
  4285.  
  4286.          whichsat=-1;
  4287.          for (i=0;i<=NumberOfSats;i++){
  4288.             nc=strlen(SATS[i].SatsName);
  4289.             if (strncmpi(SatName,SATS[i].SatsName,nc) == 0){
  4290.                whichsat=i;
  4291.             }
  4292.          }
  4293.       }
  4294.  
  4295. /*  
  4296.   if whichsat is >= 0 then check to see if we are on line
  4297.   with the word satellite (imatch=0) or on a mode data line
  4298.   (imatch=1)
  4299.  
  4300.   if satellite line (imatch=0) then get alon,alat and beacon frequency
  4301.   otherwise store mode data
  4302. */
  4303.       if(whichsat >= 0){
  4304.          if (imatch == 1){
  4305.               sscanf(line1,"%3d%3d%3s",&min,&max,&mode);
  4306.  
  4307.         Modes[whichsat][SATS[whichsat].NumModes].MinPhase = min;
  4308.             Modes[whichsat][SATS[whichsat].NumModes].MaxPhase = max;
  4309.  
  4310.         strcpy(Modes[whichsat][SATS[whichsat].NumModes].ModeStr,mode);
  4311.  
  4312.         SATS[whichsat].NumModes++;
  4313.  
  4314. /*  too many ? */
  4315.             if (SATS[whichsat].NumModes > MaxModes){
  4316.                printf (" Too many modes specified for %s\n",
  4317.                   SATS[whichsat].SatsName);
  4318.                exit(-1);
  4319.             }
  4320.      }
  4321.          if (imatch == 0){
  4322.             imatch = 1;
  4323.  
  4324.         fgets(line1,80,ModeDataFile);
  4325.         sscanf(line1,"%f%f%f",&frqbec,&Alon,&Alat);
  4326.  
  4327. /*  set squint computation option */
  4328.  
  4329.             if (Alon == ZERO && Alat == ZERO) {
  4330.                SATS[whichsat].SquintOpt=0;
  4331.             } else {
  4332.                SATS[whichsat].SquintOpt=1;
  4333.             }
  4334.  
  4335.             SATS[whichsat].Beacon=frqbec*1.0e6;
  4336.             SATS[whichsat].Alon=Alon/D2R;
  4337.             SATS[whichsat].Alat=Alat/D2R;
  4338.      }
  4339.       }
  4340.    }
  4341.    fclose (ModeDataFile);
  4342. }
  4343. void ReadNasa(void)
  4344. /*  read kepler data in nasa 2 line format */
  4345. {
  4346.    int i,nc,idum,SatNum,YearLau,LauNum;
  4347.    unsigned long int EpochRev,ElementSet;
  4348.    float EpochYear,EpochDay,DecayRate,MeanMotion;
  4349.    float Inclination,RAofAN,Eccentricity,ArgPerigee;
  4350.  
  4351. /* read until end of file */
  4352.  
  4353.    NumberOfSats = -1;
  4354.  
  4355.    while (fgets(SatName,40,NasaDataFile)){
  4356.       NumberOfSats++;
  4357.  
  4358.       if (NumberOfSats > MaxSat) {
  4359.          NumberOfSats--;
  4360.          printf (" Read in Maximum Number of Satellites %d\n",MaxSat);
  4361.          printf (" Last Satellite Read in is %s\n",
  4362.             SATS[NumberOfSats].SatsName);
  4363.          return;
  4364.       }
  4365.  
  4366. /*  find real end of this string */
  4367.       nc=RealEndOfLine(SatName);
  4368.  
  4369.       strncpy(SATS[NumberOfSats].SatsName,SatName,nc);
  4370.  
  4371. /*
  4372.       SATS[NumberOfSats].SatsName = strdup(SatName);
  4373.       SATS[NumberOfSats].SatsName[nc] = '\0';
  4374. */
  4375.  
  4376.       fgets(line1,80,NasaDataFile);
  4377.       fgets(line2,80,NasaDataFile);
  4378.  
  4379.       sscanf(line1,"%*2c%5d%*10c%2f%12f%10f%*21c%5lu",
  4380.             &SatNum,&EpochYear,&EpochDay,&DecayRate,&ElementSet);
  4381.  
  4382.       SATS[NumberOfSats].SatNum=SatNum;
  4383.  
  4384. /* strip off checksum */
  4385.       ElementSet /= 10;
  4386.       SATS[NumberOfSats].EleSet=ElementSet;
  4387.  
  4388.       sscanf(line2,"%*8c%8f%8f%7f%8f%8f%11f%5lu",
  4389.          &Inclination,&RAofAN,&Eccentricity,&ArgPerigee,&MeanAnomaly,
  4390.          &MeanMotion,&EpochRev);
  4391.  
  4392.       SATS[NumberOfSats].SquintOpt= 0;
  4393.       SATS[NumberOfSats].NumModes = 0;
  4394.  
  4395. /*  set default values for mode information */
  4396.  
  4397.       Modes[NumberOfSats][0].MinPhase = 0;
  4398.       Modes[NumberOfSats][0].MaxPhase = 256;
  4399.  
  4400.       strcpy(Modes[NumberOfSats][0].ModeStr,"---");
  4401.  
  4402. /*  set default values for beacon frequency, alat and alon */
  4403.  
  4404.       SATS[NumberOfSats].Beacon = 145.00 * 1.0e6;
  4405.       SATS[NumberOfSats].Alat   = 0.0;
  4406.       SATS[NumberOfSats].Alon   = 0.0;
  4407.  
  4408. /*  store/convert values */
  4409.  
  4410.       SATS[NumberOfSats].YE = BasicYear + EpochYear;
  4411.       SATS[NumberOfSats].TE = EpochDay;
  4412.       SATS[NumberOfSats].IN = Inclination;
  4413.       SATS[NumberOfSats].RA = RAofAN;
  4414.       SATS[NumberOfSats].EC = Eccentricity*1.0e-7;
  4415.       SATS[NumberOfSats].WP = ArgPerigee;
  4416.       SATS[NumberOfSats].MA = MeanAnomaly;
  4417.       SATS[NumberOfSats].MM = MeanMotion;
  4418.       SATS[NumberOfSats].M2 = DecayRate;
  4419.       SATS[NumberOfSats].RV = EpochRev/10;
  4420.    }
  4421. }
  4422. int RealEndOfLine(char string[])
  4423. /*  find real end of line padded with blanks */
  4424. {
  4425.    int i,last,nc;
  4426.    
  4427. /*  eat off last character...most likely carriage return */
  4428.    nc=strlen(string)-1;
  4429.  
  4430.    last=nc;
  4431.  
  4432. /*  
  4433.   start at end and work backwards...if location is blank then
  4434.   set last = i.  
  4435.  
  4436.   if string then blanks last will contain the position of the last
  4437.   character
  4438.  
  4439.   if string then blank string blank... then want to return the
  4440.   value of the first character we hit. if we have found a word
  4441.   then last contains the location of the end of the word, now if
  4442.   we hit a blank and another word last will be greater than i
  4443.   and we stop 
  4444. */
  4445.    for (i=nc;i>=0;i--){
  4446.       if(string[i] == ' '){
  4447.          last=i;
  4448.       }
  4449.       if(last > i)return last;
  4450.    }
  4451.  
  4452.    if (last == 0)last=nc;
  4453.  
  4454.    return last;
  4455. }
  4456. void ReadOtherQth(void)
  4457. /*  
  4458.     read in other qth file info.  used to do mutual visibilities 
  4459.     user can use commands (starthere) and (stophere) to pull out
  4460.     the qth(s) of interest 
  4461. */
  4462. {
  4463.    FILE *OQthFile;
  4464.  
  4465.    int On;
  4466.    char DummyLine[132];
  4467.    char *p;
  4468.  
  4469.    if ((OQthFile = fopen(DxQthFile,"r")) == 0){
  4470.       printf (" could not open file %s\n",DxQthFile);
  4471.       exit(-1);
  4472.    } 
  4473.  
  4474.    On = -1;
  4475.  
  4476.    while (fgets(DummyLine,132,OQthFile)){
  4477.  
  4478.       if(strnicmp(DummyLine,"starthere",9) == 0){
  4479.          On=1;
  4480.          fgets (DummyLine,132,OQthFile);
  4481.       }
  4482.  
  4483.       if(strnicmp(DummyLine,"stophere",8) == 0){
  4484.          On=-1;
  4485.          fgets (DummyLine,132,OQthFile);
  4486.       }
  4487.  
  4488.       if(On > 0) {
  4489.          NumberOfQth++;
  4490.  
  4491. /*  latitude */
  4492.          p=strtok(DummyLine," ");
  4493.          Qth[NumberOfQth].LA=atof(p);
  4494.  
  4495. /*  longitude */
  4496.          p=strtok(NULL," ");
  4497.          Qth[NumberOfQth].LO=atof(p);
  4498.  
  4499. /*  altitude */
  4500.          p=strtok(NULL," ");
  4501.          Qth[NumberOfQth].HT=atof(p);
  4502.  
  4503. /*  grid locator */
  4504.          p=strtok(NULL," ");
  4505.  
  4506. /*  rest of line....identifier */
  4507.          p=strtok(NULL,"\n");
  4508.  
  4509. /*  only keep first 8 characters */
  4510.          strncpy(Qth[NumberOfQth].QthLoc,p,8);
  4511.          Qth[NumberOfQth].QthLoc[8] = '\0';
  4512.       }
  4513.    } 
  4514.  
  4515.    fclose (OQthFile);
  4516.  
  4517.    printf (" Read in %d locations\n",NumberOfQth);
  4518. }
  4519. void RealTimeDisplay(int argv,char *argc[],int WhatToShow)
  4520. /*  
  4521.    real time of satellites:
  4522.       WhatToShow = 0 currently visible satellites 
  4523.       WhatToShow = 1 all satellites 
  4524.       WhatToShow = 2 selected set of satellites
  4525. */
  4526. {
  4527.    int i,nc;
  4528.    float OutputRate,ElevationLimit;
  4529.    float TimeAOS,TimeLOS,DtHour,DtMin,DtSec;
  4530.  
  4531.    OutputFile = stdout;
  4532.  
  4533. /*  
  4534.   set up array of satellites to use
  4535.  
  4536.   also catch non defined cases here.....
  4537. */
  4538.  
  4539.    switch (WhatToShow){
  4540.       case 0:
  4541.       case 1:
  4542.          NumberOfSatsSelected = NumberOfSats; 
  4543.  
  4544.          for (i=0;i<=NumberOfSats;i++){
  4545.             WhichSatelliteSelected[i] = i;
  4546.          }
  4547.  
  4548.          break;
  4549.  
  4550.       case 2:
  4551.          
  4552.          if ((SelectFile = fopen ("select.dat","r")) == 0){
  4553.             printf (" cant open select.dat\n");
  4554.             printf (" run STP with the S option to generate. (STP S)\n");
  4555.             exit(-1);
  4556.          }
  4557.  
  4558.          NumberOfSatsSelected = -1;
  4559.  
  4560.          while(fgets(SatName,40,SelectFile)){
  4561.             nc=strlen(SatName)-1;
  4562.  
  4563. /*  handle cases where null line entered....strncmpi will match all lines...*/ 
  4564.  
  4565.             if (nc > 0) {
  4566.                for (i=0;i<=NumberOfSats;i++){
  4567.                   if (strncmpi(SatName,SATS[i].SatsName,nc) == 0){
  4568.                      NumberOfSatsSelected++;
  4569.                      WhichSatelliteSelected[NumberOfSatsSelected] = i;
  4570.                   }
  4571.                }
  4572.             }
  4573.          }
  4574.  
  4575.          if (NumberOfSatsSelected <= -1){
  4576.             printf (" No sats matched from select.dat\n");
  4577.             exit(-1);
  4578.          }
  4579.  
  4580.          fclose (SelectFile);
  4581.  
  4582.          break;
  4583.  
  4584.       default:
  4585.          printf (" can't match option %d in RealTimeDisplay\n",
  4586.             WhatToShow);
  4587.          exit(-1);
  4588.     }
  4589.  
  4590. /*  initialize all values in array */
  4591.     for (i=0;i<MaxSat;i++){
  4592.        UpDown[i] = 0;
  4593.     }
  4594.  
  4595. /*  see if user specified an update rate...otherwise set default */
  4596.     OutputRate = 1.0;
  4597.  
  4598.     if (argv == 3){
  4599.        OutputRate = atof(argc[2]);
  4600.     }
  4601.  
  4602. /* set observer as local site */
  4603.    WhichObserver = 0;
  4604.  
  4605. /*  initailize Qth parameters */
  4606.  
  4607.    InitQth();
  4608.  
  4609. /*  
  4610.   initialize rise/set times
  4611.  
  4612.   performed here to minmize switching in real time loop
  4613.   also makes variable housekeeping easier
  4614. */
  4615.  
  4616.    switch (WhatToShow){
  4617.       case 0:
  4618.          break;
  4619.  
  4620.       case 1:
  4621.       case 2:
  4622.          UTCLOCTime();
  4623.  
  4624.          Year   = (float)GMTyear + BasicYear;
  4625.          Month  = (float)GMTmon + ONE;
  4626.          Day    = (float)GMTdaym;
  4627.  
  4628.          Hour   = (float)GMThour;
  4629.          Minute = (float)GMTmin;
  4630.          Second = (float)GMTsec;
  4631.  
  4632. /*  get day number */
  4633.  
  4634.          DayNow = DayNumber(Year,Month,Day);
  4635.  
  4636. /*  loop over satellites */
  4637.  
  4638.          for (i=0;i<=NumberOfSatsSelected;i++){
  4639.             WhichSatellite = WhichSatelliteSelected[i];
  4640.  
  4641. /*  current time in hours.... set here since FindAOS/FindLOS reset */
  4642.  
  4643.             TimeNow = (Hour+Minute/60.0+Second/3600.0)/24.0;
  4644.  
  4645. /*  initailize sun parameters...a must do for each satellite */
  4646.  
  4647.             InitSun();
  4648.  
  4649. /*  initailize antenna parameters...a must do for each satellite */
  4650.  
  4651.             InitAntenna();
  4652.  
  4653. /*  compute satellite position */
  4654.             SatellitePosition();
  4655.  
  4656. /* determine azimiuth/elevation between observer and satellite */
  4657.             ObserverSatelliteGeometry();
  4658. /*  
  4659.   find AOS/LOS times.... current clock time 
  4660. */
  4661.  
  4662.              if (Elevation >= ZERO){
  4663.                 RiseSet[WhichSatellite].RiseHour   = Hour;
  4664.                 RiseSet[WhichSatellite].RiseMinute = Minute;
  4665.                 RiseSet[WhichSatellite].RiseSecond = Second;
  4666.              }else{
  4667.                 TimeAOS = FindAOS();
  4668.  
  4669.                 HourMinuteSecond(TimeAOS,stdout,-999);
  4670.  
  4671.                 RiseSet[WhichSatellite].RiseHour   = TempHour;
  4672.                 RiseSet[WhichSatellite].RiseHour  += TempDay*24.0;
  4673.                 RiseSet[WhichSatellite].RiseMinute = TempMin;
  4674.                 RiseSet[WhichSatellite].RiseSecond = TempSec;
  4675.  
  4676. /* if TimeAOS is negative then rise time is greater than 2 days away */
  4677.  
  4678.                 if (TimeAOS < ZERO){
  4679.                    RiseSet[WhichSatellite].RiseHour   = Hour + 48.0;
  4680.                    RiseSet[WhichSatellite].RiseMinute = Minute;
  4681.                    RiseSet[WhichSatellite].RiseSecond = Second;
  4682.                 }
  4683.              }
  4684.  
  4685.              TimeLOS = FindLOS();
  4686.              HourMinuteSecond(TimeLOS,stdout,-999);
  4687.  
  4688.              RiseSet[WhichSatellite].SetHour   = TempHour;
  4689.              RiseSet[WhichSatellite].SetHour  += TempDay*24.0;
  4690.              RiseSet[WhichSatellite].SetMinute = TempMin;
  4691.              RiseSet[WhichSatellite].SetSecond = TempSec;
  4692.  
  4693. /* if TimeLOS is negative then set time is greater than 2 days away */
  4694.  
  4695.              if (TimeLOS < ZERO){
  4696.                 RiseSet[WhichSatellite].SetHour   = Hour + 48.0;
  4697.                 RiseSet[WhichSatellite].SetMinute = Minute;
  4698.                 RiseSet[WhichSatellite].SetSecond = Second;
  4699.              }
  4700.          }
  4701.          break;
  4702.    }
  4703.  
  4704. /*  real time loop over each satellite */
  4705.  
  4706.    while (bioskey(1) == 0){
  4707.  
  4708. /*  get current time */
  4709.       UTCLOCTime();
  4710.  
  4711.       Year   = (float)GMTyear + BasicYear;
  4712.       Month  = (float)GMTmon + ONE;
  4713.       Day    = (float)GMTdaym;
  4714.  
  4715.       Hour   = (float)GMThour;
  4716.       Minute = (float)GMTmin;
  4717.       Second = (float)GMTsec;
  4718.  
  4719. /*  get day number */
  4720.  
  4721.       DayNow = DayNumber(Year,Month,Day);
  4722.  
  4723. /*  if times up then dump data */
  4724.  
  4725.       if (DeltTime(OutputRate)){
  4726.          ClearTheScreen();
  4727.  
  4728. /*  show user current local and GMT time */
  4729.  
  4730.          printf (" LOC %10s %10s %02d 19%02d %02d:%02d:%02d\n",
  4731.             DayNames[LOCdayw],
  4732.             MonthNames[LOCmon],
  4733.             LOCdaym,
  4734.             LOCyear,
  4735.             LOChour,
  4736.             LOCmin,
  4737.             LOCsec);
  4738.  
  4739.          printf (" GMT %10s %10s %02d 19%02d %02d:%02d:%02d\n",
  4740.             DayNames[GMTdayw],
  4741.             MonthNames[GMTmon],
  4742.             GMTdaym,
  4743.             GMTyear,
  4744.             GMThour,
  4745.             GMTmin,
  4746.             GMTsec);
  4747.  
  4748.          switch (WhatToShow){
  4749.             case 0:
  4750.                printf ("  Az    El    Range Squnt  Doppler");
  4751.                printf ("  MA    Mode Sun\n");
  4752.                break;
  4753.  
  4754.             case 1:
  4755.             case 2:
  4756.                printf ("  Az    El    Range Squnt  Doppler");
  4757.                printf ("  MA    Mode Sun Status AOS--LOS\n");
  4758.                break;
  4759.          }
  4760.  
  4761.          for (i=0;i<=NumberOfSatsSelected;i++){
  4762.             WhichSatellite = WhichSatelliteSelected[i];
  4763.  
  4764. /*  current time in hours.... set here since FindAOS/FindLOS reset */
  4765.  
  4766.             TimeNow = (Hour+Minute/60.0+Second/3600.0)/24.0;
  4767.  
  4768. /*  initailize sun parameters...a must do for each satellite */
  4769.  
  4770.             InitSun();
  4771.  
  4772. /*  initailize antenna parameters...a must do for each satellite */
  4773.  
  4774.             InitAntenna();
  4775.  
  4776. /*  compute satellite position */
  4777.  
  4778.             SatellitePosition();
  4779.  
  4780. /* determine azimiuth/elevation between observer and satellite */
  4781.  
  4782.             ObserverSatelliteGeometry();
  4783.  
  4784. /*  set elevation limit based on option */
  4785.  
  4786.             switch (WhatToShow){
  4787.                case 0:
  4788.                   ElevationLimit = ZERO;
  4789.                   break;
  4790.  
  4791.                case 1:
  4792.                case 2:
  4793.                   ElevationLimit = -1000.0;
  4794.                   break;
  4795.             }
  4796.  
  4797.             if (Elevation > ElevationLimit){
  4798.  
  4799. /*  compute sun lit status */
  4800.  
  4801.                SunData();
  4802.  
  4803.                printf(" %4.0f %4.0f",
  4804.                   Azimuth*D2R,
  4805.                   Elevation*D2R);
  4806.  
  4807. /*  compute squint/rangerate */
  4808.  
  4809.                SquintRangeRate();
  4810.  
  4811.  
  4812. /* 
  4813.    which way is elevation changing....
  4814.       if RangeRate is positive then going away
  4815.       if RangeRate is negative then coming towards
  4816. */
  4817.  
  4818.                if (RangeRate >= ZERO)printf (" v");
  4819.                if (RangeRate  < ZERO)printf (" ^");
  4820.  
  4821.                printf (" %6.0f",
  4822.                   RangeToSatellite);
  4823.  
  4824. /*  get mode info */
  4825.  
  4826.                WhatMode();
  4827.  
  4828. /*  if Omnidirectional antennas on ao-10/ao-13... then adjust squint angle */
  4829.  
  4830.                if (strncmpi(CurrentMode,"Bo",2) == 0) {
  4831.                   SquintAngle = ABS(90.0 - SquintAngle);
  4832.                }
  4833.  
  4834.                printf (" %4.0f",SquintAngle);
  4835.  
  4836. /*  compute doppler */
  4837.  
  4838.                Doppler(SATS[WhichSatellite].Beacon);
  4839.  
  4840.                printf ("  %5.1f   %s",CurrentMAd,CurrentMode);
  4841.  
  4842. /*  sunlit status */
  4843.  
  4844.                printf ("%s",ECL);
  4845.  
  4846. /*  
  4847.   different options show some extra stuff
  4848.  
  4849.   also recompute AOS/LOS times if sat has just risen or set
  4850. */
  4851.  
  4852.                switch (WhatToShow){
  4853.                   case 0:
  4854.                      break;
  4855.  
  4856.                   case 1:
  4857.                   case 2:
  4858.                      if (Elevation > ZERO){
  4859. /*
  4860.                         printf ("  Up  ");
  4861. */
  4862.                         printf (" -LOS-");
  4863.  
  4864. /* if satellite was down and now up compute new LOS time */
  4865.  
  4866.                         if (UpDown[WhichSatellite] == -1){
  4867.  
  4868.                            TimeLOS = FindLOS();
  4869.                            HourMinuteSecond(TimeLOS,stdout,-999);
  4870.  
  4871.                            RiseSet[WhichSatellite].SetHour   = TempHour;
  4872.                            RiseSet[WhichSatellite].SetHour  += TempDay*24.0;
  4873.                            RiseSet[WhichSatellite].SetMinute = TempMin;
  4874.                            RiseSet[WhichSatellite].SetSecond = TempSec;
  4875.                         }
  4876.  
  4877.                         UpDown[WhichSatellite] = 1;
  4878.  
  4879.                         DtHour = RiseSet[WhichSatellite].SetHour - Hour;
  4880.                         DtMin  = RiseSet[WhichSatellite].SetMinute - Minute;
  4881.                         DtSec  = RiseSet[WhichSatellite].SetSecond - Second;
  4882.  
  4883.                         if (DtSec < ZERO){
  4884.                            DtMin -= ONE;
  4885.                            DtSec += 60.0;
  4886.                         }
  4887.                         if (DtMin < ZERO){
  4888.                             DtHour -= ONE;
  4889.                             DtMin  += 60.0;
  4890.                         }
  4891.                      }else{
  4892. /*
  4893.                         printf (" Down ");
  4894. */
  4895.                         printf ("  aos ");
  4896.  
  4897. /* if satellite was up and now down compute new AOS time */
  4898.  
  4899.                         if (UpDown[WhichSatellite] == 1){
  4900.                            UpDown[WhichSatellite] = -1;
  4901.  
  4902.                            TimeAOS = FindAOS();
  4903.                            HourMinuteSecond(TimeAOS,stdout,-999);
  4904.  
  4905.                            RiseSet[WhichSatellite].RiseHour   = TempHour;
  4906.                            RiseSet[WhichSatellite].RiseHour  += TempDay*24.0;
  4907.                            RiseSet[WhichSatellite].RiseMinute = TempMin;
  4908.                            RiseSet[WhichSatellite].RiseSecond = TempSec;
  4909.                         }
  4910.  
  4911.                         DtHour = RiseSet[WhichSatellite].RiseHour - Hour;
  4912.                         DtMin  = RiseSet[WhichSatellite].RiseMinute - Minute;
  4913.                         DtSec  = RiseSet[WhichSatellite].RiseSecond - Second;
  4914.  
  4915.                         if (DtSec < ZERO){
  4916.                            DtMin -= ONE;
  4917.                            DtSec += 60.0;
  4918.                         }
  4919.                         if (DtMin < ZERO){
  4920.                             DtHour -= ONE;
  4921.                             DtMin  += 60.0;
  4922.                         }
  4923.                      }
  4924.  
  4925. /*  AOS/LOS times.....Relative */
  4926.  
  4927.                      printf (" %02.0f:%02.0f:%02.0f",
  4928.                         DtHour,
  4929.                         DtMin,
  4930.                         DtSec);
  4931.  
  4932.                      break;
  4933.  
  4934.                }
  4935.  
  4936.                printf (" %12s\n",
  4937.                   SATS[WhichSatellite].SatsName);
  4938.             }
  4939.          }   /*  end of WhichSatellite loop */
  4940.  
  4941. /*  save current time for output determination */
  4942.  
  4943.          HourOld   = Hour;
  4944.          MinuteOld = Minute;
  4945.          SecondOld = Second;
  4946.       }
  4947.    }
  4948. }
  4949. void RunCheckSum(int argv, char *argc[])
  4950. /*  perform checksum calculations as a standalone feature */
  4951. {
  4952. /*  if user did not specify filename on command line then ask for it */
  4953.  
  4954.     if (argv <= 2){
  4955.        printf (" Filename to check: ");
  4956.        gets(line1);
  4957.     }else{
  4958.        strcpy(line1,argc[2]);
  4959.     }
  4960.  
  4961. /*  perform checksum calculations */
  4962.  
  4963.     DoCheckSum(line1,1);
  4964. }
  4965. void RevolutionCorrection(void)
  4966. /*  
  4967.    read revolution correction database...if it exists
  4968.    if it does then reset rev counter
  4969. */
  4970. {
  4971.    int i,nc;
  4972.    float RevCorrection;
  4973.  
  4974. /*  if file does not exist then return */
  4975.    if ((RevolutionDataFile = fopen(RevCorFile,"r")) == 0){
  4976.       return;
  4977.    }
  4978.  
  4979. /* 
  4980.     read file
  4981.     loop over input satelites and try to match satellite
  4982.     correct data
  4983. */
  4984.  
  4985.    while (fgets(line1,80,RevolutionDataFile)){
  4986.  
  4987.       fgets(line2,80,RevolutionDataFile);
  4988.       sscanf(line2,"%f",&RevCorrection);
  4989.  
  4990.       for (i=0;i<=NumberOfSats;i++){
  4991.          nc=strlen(line1)-1;
  4992.          if (strnicmp(line1,SATS[i].SatsName,nc) == 0){
  4993.             SATS[i].RV=SATS[i].RV + RevCorrection;
  4994.          }
  4995.       }
  4996.    }
  4997.  
  4998.    fclose (RevolutionDataFile);
  4999. }
  5000. void SatelliteDump(void)
  5001. /*  dump out info on satellites */
  5002. {
  5003.    int i;
  5004.    float CurrentDay;
  5005.  
  5006.    UTCLOCTime();
  5007.  
  5008.    Month = (float)GMTmon + ONE;
  5009.    Day   = (float)GMTdaym;
  5010.    Year  = (float)GMTyear + BasicYear;
  5011.  
  5012. /*  get current day number */
  5013.    CurrentDay=DayNumber(Year,Month,Day);
  5014.  
  5015.    printf (" #    Beacon  Alon  Alat");
  5016.    printf ("   Apogee    Perigee   Period");
  5017.    printf ("   Ecc     Age");
  5018.    printf ("\n");
  5019.  
  5020.    printf ("        Mhz   deg   deg ");
  5021.    printf ("     Km         Km       min ");
  5022.    printf ("           days");
  5023.    printf ("\n");
  5024.  
  5025.    for (i=0;i<=NumberOfSats;i++){
  5026.  
  5027.       printf (" %2d %9.4f % 4.0f % 4.0f",
  5028.       SATS[i].NumModes,
  5029.       SATS[i].Beacon*1.0e-6,
  5030.       SATS[i].Alon*D2R,
  5031.       SATS[i].Alat*D2R);
  5032.  
  5033. /*  apogee,perigee calculations */
  5034.       printf (" % 9.1f % 9.1f   % 7.1f",
  5035.          SATS[i].A0*(ONE+SATS[i].EC)-RE,
  5036.          SATS[i].A0*(ONE-SATS[i].EC)-RE,
  5037.          SATS[i].Period);
  5038.  
  5039.       printf (" % 8.5f",
  5040.          SATS[i].EC);
  5041.  
  5042.       printf (" % 6.2f",
  5043.          CurrentDay-(SATS[i].DE+SATS[i].TE));
  5044.  
  5045.       printf (" %s \n",
  5046.          SATS[i].SatsName);
  5047.    }
  5048. }
  5049. void SatellitePosition(void)
  5050. /* calculate satellite position at DayNow,TimeNow */
  5051. {
  5052.    float DT,KD,KDP,DR,EA,EccSin,EccCos,DEA,FKepler,FKeplerP;
  5053.    float CW,SW,CQ,SQ;
  5054.    float C,S,A,B;
  5055.    float Denominator;
  5056.    float tempvar;
  5057.  
  5058. /* elapsed time since epoch, days */
  5059.    TimeSinceEpoch = (DayNow  - SATS[WhichSatellite].DE) + 
  5060.                     (TimeNow - SATS[WhichSatellite].TE);       
  5061.  
  5062. /* linear drag terms */
  5063.    DT  = 0.50*SATS[WhichSatellite].DC*TimeSinceEpoch;            
  5064.    KD  = ONE + 4.0*DT;
  5065.    KDP = ONE - 7.0*DT;
  5066.  
  5067. /* mean anomaly at Year,TimeNow */
  5068.    MeanAnomaly = SATS[WhichSatellite].MA 
  5069.                 +SATS[WhichSatellite].MM*TimeSinceEpoch*(ONE-3.0*DT);  
  5070.  
  5071. /* strip out whole # of revs */
  5072.    DR = AINT(MeanAnomaly/TWOPI);          
  5073.  
  5074. /* M now in range 0 - 2 pi */
  5075.    MeanAnomaly = MeanAnomaly - DR*TWOPI;            
  5076.  
  5077. /* current orbit number */
  5078.    RevolutionNumber = SATS[WhichSatellite].RV + DR;                
  5079.  
  5080. /* 
  5081.    solve M = EA - EC*sin(EA) for EA given M by 
  5082.    Newtons Method with acceleration
  5083.  
  5084.    Inital guess for EA taken from:
  5085.  
  5086.       Fundamentals of Celestial Mechanics
  5087.       J.M.A. Danby
  5088.       2nd Edition
  5089.  
  5090. */
  5091.    EA=MeanAnomaly+0.85*SGN(sin(MeanAnomaly))*SATS[WhichSatellite].EC;
  5092.  
  5093.    EccSin=SATS[WhichSatellite].EC*sin(EA);
  5094.  
  5095.    FKepler=EA-EccSin-MeanAnomaly;
  5096.  
  5097.    while (ABS(FKepler) > 0.00001){
  5098.       EccCos   =  SATS[WhichSatellite].EC*cos(EA);
  5099.       FKeplerP =  ONE-EccCos;
  5100.       DEA      = -FKepler/FKeplerP;
  5101.       DEA      = -FKepler/(FKeplerP+0.50*DEA*EccSin);
  5102.       DEA      = -FKepler/(FKeplerP+0.50*DEA*EccSin+DEA*DEA*EccCos/6.0);
  5103.       EA      +=  DEA;
  5104.       EccSin   =  SATS[WhichSatellite].EC*sin(EA);
  5105.       FKepler  =  EA-EccSin-MeanAnomaly;
  5106.    }
  5107.  
  5108.    C = cos(EA);
  5109.    S = sin(EA);
  5110.  
  5111.    Denominator=ONE-SATS[WhichSatellite].EC*C;
  5112.    
  5113.    A = SATS[WhichSatellite].A0*KD;
  5114.    B = SATS[WhichSatellite].B0*KD;
  5115.    AltitudeOfSatellite = A*Denominator;
  5116.  
  5117. /* calculate satellite positions and velocity in plane of ellipse */
  5118.  
  5119.    tempvar=SATS[WhichSatellite].N0/Denominator;
  5120.  
  5121.    Sx = A*(C-SATS[WhichSatellite].EC);
  5122.    Vx = -A*S*tempvar;
  5123.    Sy = B*S;
  5124.    Vy = B*C*tempvar;
  5125.  
  5126.    tempvar = TimeSinceEpoch*KDP;
  5127.  
  5128.    AP = SATS[WhichSatellite].WP
  5129.        +SATS[WhichSatellite].WD*tempvar;
  5130.    CW = cos(AP);
  5131.    SW = sin(AP);
  5132.  
  5133.    RAAN = SATS[WhichSatellite].RA
  5134.          +SATS[WhichSatellite].QD*tempvar;
  5135.    CQ = cos(RAAN);
  5136.    SQ = sin(RAAN);
  5137.  
  5138. /* plane -> celestial coordinate transformation [C] = [RAAN][IN][AP] */
  5139.  
  5140.    CXx =  CW*CQ-SW*SATS[WhichSatellite].CI*SQ;
  5141.    CXy = -SW*CQ-CW*SATS[WhichSatellite].CI*SQ;
  5142.    CXz =           SATS[WhichSatellite].SI*SQ;
  5143.  
  5144.    CYx =  CW*SQ+SW*SATS[WhichSatellite].CI*CQ;
  5145.    CYy = -SW*SQ+CW*SATS[WhichSatellite].CI*CQ;
  5146.    CYz =       -   SATS[WhichSatellite].SI*CQ;
  5147.  
  5148.    CZx =        SW*SATS[WhichSatellite].SI;
  5149.    CZy =        CW*SATS[WhichSatellite].SI;
  5150.    CZz =           SATS[WhichSatellite].CI;
  5151.  
  5152. /* compute satellites position vector, antenna axis unit vector and
  5153.    velocity in celestial coordinates
  5154. */
  5155.  
  5156.    SATx = Sx*CXx+Sy*CXy;
  5157.    SATy = Sx*CYx+Sy*CYy;
  5158.    SATz = Sx*CZx+Sy*CZy;
  5159.  
  5160.    ANTx = SATS[WhichSatellite].ax*CXx
  5161.          +SATS[WhichSatellite].ay*CXy
  5162.          +SATS[WhichSatellite].az*CXz;
  5163.    ANTy = SATS[WhichSatellite].ax*CYx
  5164.          +SATS[WhichSatellite].ay*CYy
  5165.          +SATS[WhichSatellite].az*CYz;
  5166.    ANTz = SATS[WhichSatellite].ax*CZx
  5167.          +SATS[WhichSatellite].ay*CZy
  5168.          +SATS[WhichSatellite].az*CZz;
  5169.  
  5170.    VELx = Vx*CXx+Vy*CXy;
  5171.    VELy = Vx*CYx+Vy*CYy;
  5172.    VELz = Vx*CZx+Vy*CZy;
  5173.  
  5174. /* express SAT, ANT and VEL in geocentric coordinates */
  5175.  
  5176.    GHAA = GHAE + WE*TimeSinceEpoch;
  5177.  
  5178.    C = cos(-GHAA);
  5179.    S = sin(-GHAA);
  5180.  
  5181.    Sx = SATx*C - SATy*S;
  5182.    Sy = SATx*S + SATy*C;
  5183.    Sz = SATz;
  5184.  
  5185.    Ax = ANTx*C - ANTy*S;
  5186.    Ay = ANTx*S + ANTy*C;
  5187.    Az = ANTz;
  5188.  
  5189.    Vx = VELx*C - VELy*S;
  5190.    Vy = VELx*S + VELy*C;
  5191.    Vz = VELz;
  5192. }
  5193. void SatellitePredict(int argv,char *argc[],int Option)
  5194. /*  controlling routine to predict satellite passes */
  5195. {
  5196.    float CurrentDay,DayEnd,NumberOfDays;
  5197.    float HourStart;
  5198.    int i,nc;
  5199.  
  5200. /*  set observer to home Qth */
  5201.  
  5202.    WhichObserver=0;
  5203.  
  5204. /*  
  5205.   if user did not include satellite name on command line
  5206.   then show what we got and ask...
  5207. */
  5208.  
  5209.    if (argv <= 2){
  5210.       for (i=0;i<=NumberOfSats;i++){
  5211.          printf (" %s\n",SATS[i].SatsName);
  5212.       }
  5213.  
  5214.       printf (" Which Satellite from the list ?  ");
  5215.  
  5216.       gets(SatName);
  5217.    } else {
  5218.       strcpy(SatName,argc[2]);
  5219.    }
  5220.  
  5221.    nc=strlen(SatName);
  5222.  
  5223. /*  can we match it */
  5224.  
  5225.    WhichSatellite=-1;
  5226.  
  5227.    for (i=0;i<=NumberOfSats;i++){
  5228.       if (strnicmp(SatName,SATS[i].SatsName,nc) == 0){
  5229.          if (WhichSatellite == -1){
  5230.             WhichSatellite = i;
  5231.          }else{
  5232.             printf (" Multiple Names.  Please be more specific \n");
  5233.             exit(-1);
  5234.          }
  5235.       }
  5236.    }
  5237.  
  5238.    if (WhichSatellite < 0){
  5239.       printf (" Could Not Match Satellite %s in List\n",SatName);
  5240.       exit (-1);
  5241.    }else{
  5242.       printf (" Matched your input as %s\n",SATS[WhichSatellite].SatsName);
  5243.    }
  5244.  
  5245. /*  
  5246.   get local date to check input (currently only warns user of 
  5247.   back propagation) and possibly for use as input
  5248. */
  5249.    UTCLOCTime();
  5250.  
  5251.    Month = (float)GMTmon + ONE;
  5252.    Day   = (float)GMTdaym;
  5253.    Year  = (float)GMTyear + BasicYear;
  5254.  
  5255. /*  get current day number */
  5256.    CurrentDay=DayNumber(Year,Month,Day);
  5257.  
  5258.    printf (" Enter Start Date (month day year) (<cr> for current date) : ");
  5259.  
  5260. /*  allow user option to hit return to get current date */
  5261.    strcpy(line1,GetStringFromUser(&nc));
  5262.  
  5263. /*  if nc > 0 then user input...otherwise use current information */
  5264.  
  5265.    if (nc > 0){
  5266.       sscanf(line1,"%f%f%f",&Month,&Day,&Year);
  5267.    }else{
  5268.       Month = (float)GMTmon + ONE;
  5269.       Day   = (float)GMTdaym;
  5270.       Year  = (float)GMTyear + BasicYear;
  5271.    }
  5272.  
  5273. /*  check for back propagation times */
  5274.  
  5275.    if (Month < (float)LOCmon){
  5276.       printf (" ********* Input Month < Current Month\n\n");
  5277.    }
  5278.    if (Day < (float)LOCdaym && Month <= (float)LOCmon){
  5279.       printf (" ********* Input Day < Current Day\n\n");
  5280.    }
  5281.    if (Year < BasicYear +(float)LOCyear){
  5282.       printf (" ********* Input Year < Current Year\n");
  5283.       printf ("           MA calculations will be wrong\n\n");
  5284.    }
  5285.  
  5286.    printf (" Number of Days for Prediction: ");
  5287.    NumberOfDays=GetFloatNumberFromUser(&nc);
  5288.  
  5289.    printf (" Step Size (minutes....input negative value for seconds): ");
  5290.    Step=GetFloatNumberFromUser(&nc);
  5291.  
  5292. /* 
  5293.    if Step is negative then user has input time step in seconds...
  5294.    convert to minutes 
  5295. */
  5296.  
  5297.    if (Step < ZERO){
  5298.       Step = ABS(Step)/60.0;
  5299.    }
  5300.  
  5301. /*  if option P...Option = 1 then ask for start hour */
  5302.  
  5303.    if (Option){
  5304.       printf (" Start Hour (<cr> for 0.0): ");
  5305.       HourStart = GetFloatNumberFromUser(&nc);
  5306.    }else{
  5307.       HourStart = 0.0;
  5308.    }
  5309.  
  5310.    printf (" FileName to write output to (cr for screen): ");
  5311.  
  5312. /* read filename here */
  5313.  
  5314.    strcpy(FileName,GetStringFromUser(&nc));
  5315.  
  5316.    if (nc > 0){
  5317.       if ((OutputFile = fopen(FileName,"w")) == 0){
  5318.          printf ("Can't open file %s\n",FileName);
  5319.          exit(-1);
  5320.       }
  5321.    }else{
  5322.       OutputFile = stdout;
  5323.    }
  5324.  
  5325.    DayStart = DayNumber(Year,Month,Day);
  5326.    DayEnd   = DayStart + NumberOfDays - ONE;
  5327.  
  5328. /*  initailize Qth parameters */
  5329.  
  5330.    InitQth();
  5331.  
  5332. /*  initailize sun parameters */
  5333.  
  5334.    InitSun();
  5335.  
  5336. /*  initailize antenna parameters */
  5337.  
  5338.    InitAntenna();
  5339.  
  5340. /*  footprint calculations.... */
  5341.  
  5342.    if (PlotFlag){
  5343.       FootPrintInit();
  5344.    }
  5345.  
  5346. /*  loop over time */
  5347.  
  5348.    RunningCounter    = 0;
  5349.    RunningCounterOld = 0;
  5350.  
  5351.    for (DayNow=DayStart;DayNow<=DayEnd;DayNow++){
  5352.  
  5353.       if (DayNow > DayStart)HourStart = ZERO;
  5354.  
  5355.       for (Hour=HourStart;Hour<=23.0;Hour++){
  5356.          for (Minute=ZERO;Minute<=60.0-Step;Minute=Minute+Step){
  5357.  
  5358.             RunningCounter++;
  5359.  
  5360.             TimeNow = (Hour+Minute/60.0)/24.0;
  5361.  
  5362. /*  get satellite location */
  5363.             SatellitePosition();
  5364.  
  5365. /*  determine azimuth and elevation from observers location */
  5366.             ObserverSatelliteGeometry();
  5367.  
  5368. /*  
  5369.    can we see it....
  5370.    if we can compute sun illumination 
  5371.    squint angle (antenna pointing)
  5372.    ..... print data
  5373. */
  5374.             if (Elevation >= LOSAngle(Azimuth)){
  5375.                SunData();
  5376.                SquintRangeRate();
  5377.                DataPrint(CurrentDay);
  5378.  
  5379.                if (PlotFlag){
  5380.                   FootPrint();
  5381.                }
  5382.             }
  5383.          }
  5384.       }
  5385.    }
  5386.  
  5387.    fclose (OutputFile);
  5388. }
  5389. void SelectSatellites(void)
  5390. /*  select satellites for option s */
  5391. {
  5392.    int i,ndigit;
  5393.  
  5394.    NumberOfSatsSelected=-1;
  5395.  
  5396.    if((SelectFile=fopen("select.dat","w")) == 0){
  5397.       printf (" cant open select.dat\n");
  5398.       exit(-1);
  5399.    }
  5400.  
  5401.    for (i=0;i<=NumberOfSats;i++){
  5402.       printf (" %2d %s\n",i,SATS[i].SatsName);
  5403.    }
  5404.  
  5405.    ndigit=1;
  5406.    while (ndigit > 0){
  5407.       printf (" Which Satellite from the list (<cr> to stop) ?  ");
  5408.  
  5409.       i=GetIntNumberFromUser(&ndigit);
  5410.  
  5411.       if (ndigit > 0){
  5412.          NumberOfSatsSelected++;
  5413.          WhichSatelliteSelected[NumberOfSatsSelected]=i;
  5414.          fprintf (SelectFile,"%s\n",SATS[i].SatsName);
  5415.       }
  5416.    }
  5417.  
  5418.    fclose (SelectFile);
  5419. }
  5420. void SquintRangeRate(void)
  5421. /*  compute squint angle and range rate..used for doppler */
  5422. {
  5423. /*
  5424.    resolve antenna vector along unit range vector -r.a = cos(sq)
  5425.    (SATS[].SquintOpt = 1)...e.g. AO-13, AO-10, Arsene,...
  5426.  
  5427.    assume antenna axis is aligned with satellite position vector
  5428.    (SATS[].SquintOpt = 0)...e.g. RS-10/11, AO-21, FO-20,...
  5429. */
  5430.  
  5431.    if (SATS[WhichSatellite].SquintOpt == 1) {
  5432.       SquintAngle = acos(-(Ax*Rx + Ay*Ry + Az*Rz)   )*D2R;
  5433.    } else {
  5434.       SquintAngle = acos( (Sx*Rx + Sy*Ry + Sz*Rz)/AltitudeOfSatellite)*D2R;
  5435.    }
  5436.  
  5437. /* 
  5438.   resolve sat-obs velocity vector along unit range vector
  5439.   
  5440.   note:   Qth[WhichObserver].V0z == ZERO
  5441. */
  5442.  
  5443.    RangeRate = (Vx-Qth[WhichObserver].V0x)*Rx
  5444.               +(Vy-Qth[WhichObserver].V0y)*Ry
  5445.               +(Vz                       )*Rz;
  5446. }
  5447. void SunAzimuthElevation(int PrintOut)
  5448. /*  compute sun azimuth and elevation from observers location */
  5449. {
  5450.    float MAS,TAS,C,S,Hx,Hy,Hz,SNx,SNy,SNz;
  5451.    float xRx,xRy,xRz,xR,xU,xE,xN;
  5452.  
  5453. /*  ma of sun round its orbit */
  5454.    MAS = MASE + (MASD*TimeSinceEpoch)/D2R; 
  5455.  
  5456. /*  suns true anomaly */
  5457.    TAS = MRSE + WW*TimeSinceEpoch
  5458.         +EQC1*sin(MAS)
  5459.         +EQC2*sin(2.0*MAS)
  5460.         +EQC3*sin(3.0*MAS);
  5461.  
  5462.    C = cos(TAS);
  5463.    S = sin(TAS);
  5464.  
  5465.    DistanceToSun = AstroUnit*(1.00014-0.01671*C-0.00140*cos(2.0*TAS));
  5466.  
  5467. /*  sun unit vector celestial coords */
  5468.    SNx = (float)DistanceToSun*C;
  5469.    SNy = (float)DistanceToSun*S*CNS;
  5470.    SNz = (float)DistanceToSun*S*SNS;
  5471.  
  5472.    C = cos(-GHAA);
  5473.    S = sin(-GHAA);
  5474.  
  5475. /*  sun unit vector in geocentric coordinates */
  5476.    Hx = SNx*C - SNy*S;
  5477.    Hy = SNx*S + SNy*C;
  5478.    Hz = SNz;
  5479.  
  5480. /*  azimuth/elevation */
  5481.  
  5482.    xRx = Hx-Qth[WhichObserver].Ox;
  5483.    xRy = Hy-Qth[WhichObserver].Oy;
  5484.    xRz = Hz-Qth[WhichObserver].Oz;
  5485.    xR = sqrt(xRx*xRx+xRy*xRy+xRz*xRz);
  5486.  
  5487.    xRx = xRx/xR;
  5488.    xRy = xRy/xR;
  5489.    xRz = xRz/xR;
  5490.  
  5491.    xU = xRx*Qth[WhichObserver].Ux
  5492.        +xRy*Qth[WhichObserver].Uy
  5493.        +xRz*Qth[WhichObserver].Uz;
  5494.    xE = xRx*Qth[WhichObserver].Ex
  5495.        +xRy*Qth[WhichObserver].Ey;
  5496.    xN = xRx*Qth[WhichObserver].Nx
  5497.        +xRy*Qth[WhichObserver].Ny
  5498.        +xRz*Qth[WhichObserver].Nz;
  5499.  
  5500.    CurrentSunAz = ATAN2(xE,xN)*D2R;
  5501.    CurrentSunEl = asin(xU)*D2R;
  5502.  
  5503.    if (PrintOut){
  5504.       printf (" Sun : %7.2f %7.2f %15.3f",
  5505.          CurrentSunAz,
  5506.          CurrentSunEl,
  5507.          DistanceToSun);
  5508.    }
  5509. }
  5510. void SunData(void)
  5511. /*  compute sun data */
  5512. {
  5513.    float CUA,UMD,MAS,TAS,C,S,sunt1;
  5514.    float Hx,Hy,Hz,Distance;
  5515.    float booger;
  5516.  
  5517. /*  ma of sun round its orbit */
  5518.    MAS = MASE + (MASD*TimeSinceEpoch)/D2R; 
  5519.  
  5520. /*  suns true anomaly */
  5521.    TAS = MRSE + WW*TimeSinceEpoch
  5522.         +EQC1*sin(MAS)
  5523.         +EQC2*sin(2.0*MAS)
  5524.         +EQC3*sin(3.0*MAS);
  5525.  
  5526.    C = cos(TAS);
  5527.    S = sin(TAS);
  5528.  
  5529.    Distance = 1.00014-0.01671*cos(MAS)-0.00140*cos(2.0*MAS);
  5530.  
  5531. /*  sun unit vector celestial coords */
  5532.    SUNx = Distance*C;
  5533.    SUNy = Distance*S*CNS;
  5534.    SUNz = Distance*S*SNS;
  5535.  
  5536. /*  sin of sun angle */
  5537.    SinOfSunAngle = -(ANTx*SUNx + ANTy*SUNy + ANTz*SUNz);
  5538.  
  5539. /*  illumination */
  5540.  
  5541.    booger=ONE-SinOfSunAngle*SinOfSunAngle;
  5542.    if(booger < ZERO)booger = ZERO;
  5543.  
  5544.    Illumination = 100.0*sqrt(booger);
  5545.  
  5546. /*  cos of umbral angle */
  5547.    CUA = -(SATx*SUNx + SATy*SUNy + SATz*SUNz)/AltitudeOfSatellite;
  5548.  
  5549. /*  umbral distance in units of earth radii */
  5550.    booger=ONE-CUA*CUA;
  5551.  
  5552.    if (booger < ZERO)booger = ZERO;
  5553.  
  5554.    UMD = AltitudeOfSatellite*sqrt(booger)/RE;
  5555.  
  5556. /* 
  5557.     + means for shadow side
  5558.     - means for sunny side
  5559. */
  5560.    if (CUA >= ZERO){
  5561.       strcpy(ECL,"  +  ");
  5562.    }else{
  5563.       strcpy(ECL,"  -  ");
  5564.    }
  5565.  
  5566.    if (UMD <= ONE && CUA >= ZERO) strcpy(ECL," ECL ");
  5567.  
  5568.    C = cos(-GHAA);
  5569.    S = sin(-GHAA);
  5570.  
  5571. /*  sun unit vector in geocentric coordinates */
  5572.    Hx = SUNx*C - SUNy*S;
  5573.    Hy = SUNx*S + SUNy*C;
  5574.    Hz = SUNz;
  5575.  
  5576.    sunt1 = Hx*Qth[WhichObserver].Ux
  5577.           +Hy*Qth[WhichObserver].Uy
  5578.           +Hz*Qth[WhichObserver].Uz;
  5579.  
  5580. /*  
  5581.   if sun more than 10 degrees below horizon then satellite possibly
  5582.   is visible.  depends upon satellite signature and brightness of
  5583.   background.  (10 degrees ~ 0.17 radians)
  5584. */
  5585.  
  5586.    if (sunt1 <= -0.17 && (strcmpi(ECL," ECL ") != 0)){
  5587.       strcpy(ECL," vis ");
  5588.    }
  5589.  
  5590. /*  obtain sun unit vector in orbit coordinates */
  5591.    Hx = SUNx*CXx + SUNy*CYx + SUNz*CZx;
  5592.    Hy = SUNx*CXy + SUNy*CYy + SUNz*CZy;
  5593.    Hz = SUNx*CXz + SUNy*CYz + SUNz*CZz;
  5594.  
  5595. /*  cluge here for now to limit values */
  5596.    if (Hz > ONE)Hz = ONE;
  5597.    if (Hz < -ONE)Hz = -ONE;
  5598.  
  5599.    SunElevation = asin(Hz);
  5600.    SunAzimuth   = ATAN2(Hy,Hx);
  5601. }
  5602. float TemperatureOfSky(float Frequency)
  5603. /* 
  5604.    approximation to sky temperature as a function of frequency 
  5605.  
  5606.    Frequency (Mhz)
  5607. */
  5608. {
  5609.       return exp(10.97 - 1.53*log(Frequency));
  5610. }
  5611. void UTCLOCTime(void)
  5612. /*  
  5613.    get local time (24 hour clock) and UTC 
  5614.  
  5615.    assumes variable TZ has been set
  5616. */
  5617. {
  5618.  
  5619.    time_t t;
  5620.    struct tm *gmt,*loc;
  5621.  
  5622.    t = time(NULL);
  5623.  
  5624.    loc = localtime(&t);
  5625. /*
  5626.    printf("Local time is: %s", asctime(loc));
  5627. */
  5628.  
  5629.    LOCmon =loc->tm_mon;
  5630.    LOCdaym=loc->tm_mday;
  5631.    LOCdayw=loc->tm_wday;
  5632.    LOCdayy=loc->tm_yday;
  5633.    LOCyear=loc->tm_year;
  5634.    LOChour=loc->tm_hour;
  5635.    LOCmin =loc->tm_min;
  5636.    LOCsec =loc->tm_sec;
  5637.  
  5638.    gmt = gmtime(&t);
  5639. /*
  5640.    printf("GMT is:        %s", asctime(gmt));
  5641. */
  5642.  
  5643.    GMTmon =gmt->tm_mon;
  5644.    GMTdaym=gmt->tm_mday;
  5645.    GMTdayw=gmt->tm_wday;
  5646.    GMTdayy=loc->tm_yday;
  5647.    GMTyear=gmt->tm_year;
  5648.    GMThour=gmt->tm_hour;
  5649.    GMTmin =gmt->tm_min;
  5650.    GMTsec =gmt->tm_sec;
  5651. }
  5652. void WhatMode(void)
  5653. /* determine current mode and convert Mean Anomaly to 0-256 */
  5654. {
  5655.    int i;
  5656.  
  5657.    CurrentMAd = MeanAnomaly*128.0/PI;
  5658.    CurrentMA = AINT(CurrentMAd);
  5659.  
  5660. /*  handle back propagation cases.... */
  5661.  
  5662.    if (CurrentMAd < ZERO){
  5663.       strcpy(CurrentMode,Modes[WhichSatellite][0].ModeStr);
  5664.       return;
  5665.    }
  5666.  
  5667.    for (i=0;i<=SATS[WhichSatellite].NumModes;i++){
  5668.       if (CurrentMA >= Modes[WhichSatellite][i].MinPhase && 
  5669.           CurrentMA < Modes[WhichSatellite][i].MaxPhase)
  5670.          strcpy(CurrentMode,Modes[WhichSatellite][i].ModeStr);
  5671.    }
  5672. }
  5673.